AIMETA 2005 


Firenze, 11-15 Settembre 2005 


Volume Il 


Firenze 
University 


ATTI 
- 19 - 


y 
Ja nego 


AIMETA 2005 


Atti del XVII Congresso dell'Associazione Italiana 
di Meccanica Teorica e Applicata 


Firenze, 11-15 settembre 2005 


a cura di 


Claudio Borri 
Luca Facchini 
Giorgio Federici 
Mario Primicerio 


Volume Il 


Firenze University Press 
2006 


AIMETA 2005 : atti del XVII Congresso dell’Associazione italiana di 
meccanica teorica e applicata : Firenze, 11-15 settembre 2005 : volume II 
/ a cura di Claudio Borri, Luca Facchini, Giorgio Federici, Mario Primicerio. 
- Firenze, Firenze university press, 2006 

(Atti, 19) 

http://digital.casalini.it/8884534593 

Stampa a richiesta disponibile su http://epress.unifi.it 


ISBN 10: 88-8453-459-3 (online) 
ISBN 13: 978-88-8453-459-0 (online) 
ISBN 10: 88-8453-460-7 (print) 
ISBN 13: 978-88-8453-460-6 (print) 


531 (ed. 20) 


Meccanica-Congressi-Firenze-2005 


© 2006 Firenze University Press 


Università degli Studi di Firenze 
Firenze University Press 

Borgo Albizi, 28, 50122 Firenze, Italy 
http://epress.unifi.it/ 


Printed in Italy 


XVII Congresso AIMETA di Meccanica Teorica e Applicata 


Comitato Scientifico: 


Gianni Bartoli (Università degli Studi di Firenze) 

Davide Bigoni (Università degli Studi di Trento) 

Guido Borino (Università degli Studi di Palermo) 

Ennio Carnevale (Università degli Studi di Firenze) 

Alberto Corigliano (Politecnico di Milano) 

Massimiliano Lucchesi (Università degli Studi di Firenze) 
Paolo Luchini (Università degli Studi di Salerno) 

Aleramo Lucifredi (Università degli Studi di Genova) 

Angelo Luongo (Universita degli Studi di l Aquila) 

Ettore Pennestrì (Università degli Studi di Roma ‘Tor Vergata’) 
Mario Primicerio, Presidente (Università degli Studi di Firenze) 
Terenziano Raparelli (Università degli Studi di 1’ Aquila) 

Paolo Rissone (Università degli Studi di Firenze) 

Giampiero Spiga (Università degli Studi di Parma) 


Comitato Organizzatore: 


Franco Angotti 

Ignazio Becchi 

Claudio Borri, Presidente 
Silvia Briccoli Bati 

Carlo Cinquini, Segretario AIMETA 
Paolo Citti 

Luca Facchini, Segretario 
Giorgio Federici, Tesoriere 
Giovanni Frosali 
Francesco Martelli 

Paolo Toni 

Giovanni Vannucchi 
Andrea Vignoli 


Segreteria del Congresso: 


Gabriella Montagnani 

Chiara Serpieri 

Veronika Sustik 

Ufficio Relazioni Esterne 
Facoltà di Ingegneria 

Università degli Studi di Firenze 
Via di S. Marta, 3 

50139 Firenze 

Tel: +39 055 4796491 

Fax: +39 055 4796544 

E-mail: aimeta2005@ing.unifi.it 


L’AIMETA, Associazione Italiana di Meccanica Teorica e Applicata, costituita nel 1966, riunisce i cultori 
della Meccanica nei suoi vari indirizzi: Meccanica Generale, Meccanica dei fluidi, Meccanica delle macchine, 
Meccanica dei solidi e Meccanica delle strutture. 

Attraverso Congressi e incontri e con la rivista Meccanica, 1” Associazione si propone di stabilire contatti 
fra ricercatori che operano nei diversi indirizzi, favorendo la collaborazione ed il confronto fra conoscenze 
ed esperienze diverse. 


Consiglio Direttivo dell’ AIMETA: 


Giuliano Augusti (Vice-presidente), Roberto Bassani, Gianfranco Capriz (Past-President), Carlo Cinquini 
(Segretario), Mario di Paola, Angelo Morro (Presidente), Maurizio Pandolfi (Tesoriere). 


Il XVII Congresso AIMETA di Meccanica Teorica e Applicata si è svolto con il patrocinio di: 


AIMETA - Associazione Italiana di Meccanica Teorica e Applicata 

Comune di Firenze 

CRIACIV - Centro di Ricerca Interuniversitario di Aerodinamica delle Costruzioni e Ingegneria del vento 
Facoltà di Ingegneria di Firenze - Dipartimento di Ingegneria Civile 

Istituto e Museo di Storia della Scienza 

Università degli Studi di Firenze 


ed è stato realizzato grazie al contributo di: 


qa E 
D v CASSA IM RISPARMES 
DI FIRENZE 


AT 


JT 


SAVITRANSPORT 


Universita degli Studi di Firenze 


Indice 


Nota introduttiva di Claudio Borri 


Frattura interlaminare secondo il modo I in un laminato composito 
Stefano Bennati, Massimiliano Colleluori, Domenico Corigliano, 
Paolo Sebastiano Valvo 


Dinamica del vitreo oculare indotta dai movimenti saccadici 
Chiara Cafferata, Rodolfo Repetto, Alessandro Stocchino 


Simulation of the three-dimensional flow around a square cylinder between parallel 
walls at moderate Reynolds numbers 
Simone Camarri, Maria Vittoria Salvetti, Marcelo Buffoni, Angelo Iollo 


Rans Solutions for the Numerical Prediction of Separated Flows 
Carlo de Nicola, Benedetto Mele, Renato Tognaccini 


Utilizzo di tecniche possibilistiche nella meccanica delle strutture 
Stefano Gabriele, Claudio Valente, Fabio Brancaleoni 


Instability Characteristics of Harmonic Disturbances in a Turbulent Separation Bubble 
Astrid H. Herbst, Steve Deubelbeiss, Saskia Speer, Ardeshir Hanifi, 
Dan S. Henningson 


U-RANS Simulations Around Bluff Bodies 
Claudio Marongiu, Pier Luigi Vitagliano, Francesco Capizzano, Pietro Catalano 


Analisi dinamica deterministica ed aleatoria di oscillatori che percorrono travi 
su suolo viscoelastico 


Giuseppe Muscolino, Alessandro Palmeri 


Indici 


IX 


XI 


13 


23 


35 


45 


57 


67 


79 


91 


Nota introduttiva 


Questo post-scriptum agli atti del XVII Congresso AIMETA, raccoglie alcuni contributi, che, 
pur analizzati ed approvati dal Comitato Scientifico, non hanno potuto trovare posto nel volume 
pubblicato a settembre 2005, per motivi esclusivamente editoriali e temporali. Ciò nonostante, il 
Comitato Organizzatore, avendo preso un preciso impegno con gli Autori, ha tenuto fede a tale 
impegno, pubblicando questo addendum. Forse gli Autori avrebbero sperato in tempi più stretti, ma 
anche in questa occasione i manoscritti sono stati consegnati con un qualche ritardo. 


Intendo confermare qui quanto espresso nella mia introduzione al volume di settembre 2005 e 
cioè: “Il presente volume riunisce l’impressionante quantità di contributi selezionati e raggruppati 
nei vari settori classici della meccanica teorica ed applicata, e provenienti da una vastissima comu- 
nità scientifica: a tutti gli Autori va il sincero ringraziamento del Comitato Organizzatore del Con- 
vegno. A tali settori classici si sono aggiunti temi di valenza interdisciplinare di notevole interesse e 
di alto contenuto innovativo: per questi sono stati proposti dei Minisimposi organizzati e coordinati 
dai promotori, ai quali desidero esprimere un particolare e sentito grazie.” Adesso vorrei esplicita- 
mente ringraziare gli Autori dei lavori pubblicati in questo addendum, per il loro interesse e la loro 
collaborazione. 


Anche a nome dell’intero Comitato Organizzatore, desidero rinnovare uno specifico riconosci- 
mento ai colleghi Prof. L. Facchini (Segretario del C.O.) ed alla Sig.ra G. Montagnani, che hanno 
sostenuto il maggior carico di lavoro dell’organizzazione della stampa dell’addendum. 

Infine, ritengo doveroso ringraziare il nostro “Publisher”, la Firenze University Press, in partico- 


lare la D.ssa C. Bullo, per la comprensione, il sostegno e l’aiuto qualificato a nostro lavoro. 


Claudio Borri, Prof. Ing., Dr.-Ing. h.c. 
Presidente del Comitato Organizzatore 
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Parole chiave: materiali compositi, delaminazione, frattura, interfaccia elastica, provino DCB 


SUMMARY: Si propone un modello meccanico per il provino DCB (Double Cantilever Beam), 
comunemente utilizzato per la determinazione sperimentale della resistenza alla frattura interlaminare 
secondo il modo I nei laminati compositi. A tal fine, si schematizza il provino come l’assemblaggio di 
due sublaminati collegati fra loro da un’interfaccia elasto-fragile. Grazie alla relativa semplicità del 
modello, è possibile ricavare esplicitamente gli spostamenti dei due sublaminati, le tensioni 
interlaminari e, conseguentemente, la velocità di rilascio dell’energia potenziale totale del sistema. 
L’adozione di un opportuno criterio di frattura consente, quindi, di valutare il valore critico del carico o 
dello spostamento di estremità, come funzioni esplicite dell’ampiezza della zona delaminata. Il 


confronto con alcuni risultati sperimentali presenti in letteratura appare molto buono. 


1. INTRODUZIONE 

Un ostacolo severo alla crescente diffusione dei materiali compositi nelle diverse applicazioni 
dell’ingegneria strutturale è costituito dalla loro elevata sensibilità alla presenza di difetti ed ai 
fenomeni di degrado. La delaminazione, ovvero la frattura interlaminare tra gli strati di un laminato 
composito, è una delle modalità di crisi più comuni e più insidiose. Per questo, numerosi studi sono 
stati dedicati sia alla valutazione sperimentale della resistenza a frattura dei laminati sia alla 
modellazione dei diversi aspetti coinvolti nel fenomeno (Garg [1988]). 

Nel caso della frattura secondo il modo I, la resistenza a frattura è comunemente valutata mediante 
prove di carico su particolari provini detti “DCB” (Double Cantilever Beam), i cui risultati sono 
convenzionalmente interpretati sulla base di un modello elementare che vede il laminato come l’unione 


di due travi a mensola (fig. 1). 
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Figura 1 — Provino DCB. 


Il modello meccanico che qui si propone rappresenta un’estensione di tale modello elementare, 
dove il provino DCB è schematizzato come l’assemblaggio di due sublaminati elastici collegati fra loro 
da un’interfaccia. Quest ultima è a sua volta modellata come una distribuzione continua ed uniforme di 
molle elasto-fragili (Allix e Ladevèze [1992]). 

Il modello è sufficientemente semplice da consentire la soluzione esplicita del problema, attraverso 
il metodo delle condizioni iniziali (Hetényi [1946]). Sulla base delle sole caratteristiche meccaniche 
della matrice e delle fibre del laminato, si possono allora valutare gli spostamenti dei due sublaminati, 
le tensioni interlaminari e, conseguentemente, la velocità di rilascio dell’energia potenziale totale del 
sistema al crescere dell’ampiezza della regione delaminata. L’adozione di un opportuno criterio di 
frattura consente quindi di ricavare il valore critico del carico trasversale e dello spostamento 
dell’estremità libera del provino come funzioni esplicite dell’ampiezza della zona delaminata. Tra 
l’altro, ciò consente, a differenza di altri modelli presenti in letteratura, di valutare il carico di “prima 
delaminazione”, ovvero l’intensità del carico capace di generare la nucleazione della frattura in un 
provino inizialmente integro (Hwu et al. [1995], Reedy et al. [1997]). 

Le previsioni teoriche del modello sono state confrontate con alcuni risultati sperimentali presenti 
in letteratura (Laksimi et al. [1991], Zou et al. [2003]). Considerata la semplicità del modello proposto, 
l’accordo ottenuto appare sorprendente. 


2. UN MODELLO PER LA DELAMINAZIONE SECONDO IL MODO I 

2.1. Il problema d’equilibrio 
Il modello che si propone prevede di rappresentare il provino DCB (o meglio, una sua metà, 
considerata la simmetria del problema) come una trave elastica su supporto elastico incastrata ad un 
estremo e soggetta ad un carico P all’altro estremo (fig. 2). La trave ha lunghezza complessiva /, 
mentre la porzione delaminata, priva del supporto elastico, è lunga a. È conveniente introdurre il 
rapporto adimensionale œ= a / I, cosicché le due porzioni del provino, integra e delaminata, avranno 
rispettivamente lunghezze (1 — @ l e ad. L'assenza di delaminazione corrisponde ad a = 0. 

Il valore della costante delle molle distribuite che schematizzano il supporto elastico, k, è ricavato 
direttamente dal modulo di Young della matrice, E,,, tramite la relazione k = E,, B / H, dove B è la 
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larghezza del provino ed H è il suo semispessore. In assenza di dati sulla matrice, è ragionevole porre 
k = E, B / H, dove E; è il modulo di Young della lamina in direzione trasversale. 


Figura 2 — Il modello proposto. 


Il problema d’equilibrio è scomposto in due sottoproblemi corrispondenti ai tratti AB e BC (fig. 3). 
Inoltre, nell’ipotesi di comportamento lineare, si può applicare il principio di sovrapposizione degli 
effetti e risolvere il problema per il tratto AB considerando separatamente il carico P e la coppia Pal. 


(1 -a)l | | al 


Figura 3 — Scomposizione del problema d’equilibrio. 


Siano E = E, il modulo di Young nella direzione longitudinale e J = BH’ / 12 il momento d'inerzia 
del semilaminato considerato nel modello. L'equazione della linea elastica per il tratto AB è allora 


Elv',(s)+kv(s)=0, sel0, (1-01), (1) 
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che ammette la soluzione generale nella forma 
Vig(5) =e" Ic, cos(us) + C, sin(ys)|+ er Ic, cos(us) +C, sin(ys)], (2) 


dove si è posto 


[k 
PON EJ O) 


e Cy, Co, C3 e C4 sono costanti d'integrazione. Secondo il metodo delle condizioni iniziali (Hetényi 
[1946]), differenziando la (2), si ricavano lo spostamento, la rotazione, il momento flettente e la forza 
di taglio nel punto A(s=0): 


Vo =V 14 (0) =C, +Cy; 


Ù = V 4p (0) = WC, +C, - C, +C,); 
(4) 
M, =-EJv 4,(0) = 24? EJ (-C, +C,4); 


T, =—EJv ,,(0) =2p*EJ(C, -C,-C,¿-C,4). 


Ricavando le costanti d'integrazione dalle (4) ed utilizzando le note relazioni tra funzioni 
esponenziali e funzioni iperboliche, si può porre l’espressione della linea elastica (2) nella forma 


1 1 1 
Vig (9) =V,F, (us) + —0,F, (us) - —— MF; (us) - ——T, F, (hs) , 5 
AB of Ht po WEI 0434 vo (5) 


dove si sono introdotte le funzioni ausiliarie: 


F, (us) = cosh(yus) cos(us); 
F, (us) = > [coshiys) sin(ys) + sinh(ys) cos(us)); 


1, (6) 
F, (us) = give) sin(us); 


F,(us)= lcosh(us) sin(us) — sinh(y 5) cos(ps)]. 
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S’introducano le condizioni al bordo agli estremi A e B: 


v(0)=0; — EJv"[(1- œ)l] = œP; 
A) B) (Ta, b) 
v'(0)=0; — EJv""[A— 001] =-P. 


Tramite le (7a) e le (4), l’espressione (5) si semplifica nella seguente 


1 1 
væl) == eg Mo o FS) » (8) 


la quale, differenziata e sostituita nelle (7b), permette di ricavare le caratteristiche di sollecitazione 
all’estremo A: 


M. = -P/ u _tanh[( — 0) 1] + tan[d - o) al] + 2ogul | (9a) 
°  tanh?[(1— a) ul] — tan? [A — o) ul] — 2 cosh[(1- 0) 1] cos[(1 — a) ul] i 
fe 2P ftanh[(- oul] — tani (l - a) ¿al ]jogul + 1 Ob) 
° tanh’[(1- @) ul] - tan?[(1 — o) ul] 2 cosh[(1 — œ)ul]cos[( — a) ul] i 
Considerando l’equilibrio del tratto di trave BC, non vincolato dal letto di molle, si ricava: 
a, se [0, al], (10) 


6EJ 2EJ 


dove C; e Cs sono due costanti d'integrazione, le cui espressioni si ottengono imponendo la continuità 
dello spostamento e della rotazione nella sezione B: 


C, = Fi fcosht( — 00) al sin[(1 — af) ul] + sinh[(1 — @) 1] cos[(1 — @) ul] + 


% (11a) 
+—sinh[1-@2]sin[1-@]: 
24° 
Mo . 
C, = —simh[(l — 0) ]sin[(1— a) ul] + 
24 
5 (11b) 
a o {cosh[(1 — 0) ]sin[(1— a) gl] —simh[(1 — @) 41] cosi- æu}. 
a 
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Infine, si può ricavare lo spostamento trasversale della sezione di applicazione del carico, 


P 343 1 
== al — C.al+C,). 12 
ve 355 Ej Ato a) 


2.2. Aspetti energetici 
L’energia potenziale totale del sistema, 


Nn=U-L,, (13) 


e 


è definita come la differenza tra l'energia di deformazione elastica immagazzinata dall’interfaccia 
elastica e dalla struttura, 


(14) 


PA E © e 1 M?, (5)? Pa E Mjc (sy)? Pi 
2 EJ , 2 El 02 El 


ed il lavoro delle forze esterne, 
L= Pvg. (15) 


Sostituendo nella (14) le espressioni trovate in precedenza e sviluppando gli integrali indicati, 
sarebbe possibile calcolare esplicitamente l’energia. Tuttavia, per evitare lunghi e laboriosi sviluppi 
matematici, è conveniente sfruttare il Teorema di Clapeyron, secondo cui il lavoro virtuale delle forze 
esterne è uguale al doppio dell’energia di deformazione elastica, L, = 2U. L’energia potenziale totale 
del sistema allora diventa 


1 1 1 
II=-L-L,=--L=--Pv.. 16 
2 e e 2 e 2 E. ( ) 
Si definisce, inoltre, la velocità di rilascio dell’energia, 
G=-—=--—-, (17) 


dove a è la lunghezza della parte delaminata. Sostituendo la (12) nella (16) e questa nella (17), dopo 
alcuni passaggi qui omessi per brevità, si ricava l’espressione analitica della velocità di rilascio, 


PRE si 1 E. a a a l (18) 
2EJu 2 cosld -œu +coshl(1— a) zal} 
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È facile verificare che facendo il limite per k tendente all’infinito, ovvero considerando rigida la 
parte di provino non delaminata, l’espressione (18) s'identifica con quella riportata in letteratura per il 
semplice modello a due mensole, di cui il modello proposto può essere considerato un’estensione. 


2.3. Crescita della delaminazione 
Il criterio di crescita comunemente applicato in questo tipo di problemi afferma che la delaminazione si 
propaga quando la velocità di rilascio dell’energia attinge un valore critico Gcr. Poiché nel modello, 
grazie alla simmetria del problema, si è considerata metà struttura, è necessario moltiplicare per due e 
dividere per la larghezza del provino i valori forniti dalla (18), prima di poterli confrontare con i valori 
riportati in letteratura. Pertanto, il criterio di crescita diventa: 


Gror =— = Ger (19) 


Si può determinare, allora, il valore del carico corrispondente alla propagazione: 


HYG pBEJ 
y + 1 sin[2(— a) Ml] + sinh[2(1- ul] 
2 cos[(1—0r) ul]? + cosh[d — ar) ul Y 


p= (20) 


3. APPLICAZIONE 
A titolo di esempio, si considera il provino descritto nel lavoro di Laksimi et al. [1991] (carbonio- 
epossidica con fibre di carbonio Toray T300 e resina VICOTEX M10) avente le seguenti caratteristiche 
geometriche e meccaniche: 


l= 100 mm, B=20 mm, H=6.6 mm, 
E; = 123680 N/mm’, E,= 17978.4 N/mm’, Gcr = 0.419 N/mm. 


Il valore di k si assume pari al rapporto fra il modulo di Young E, del laminato e la lunghezza 
caratteristica nella stessa direzione, cioè il semispessore del provino, k = B E, / H = 54480 N/mm”. Le 
figure 4a e 4b mostrano l’andamento della velocità di rilascio dell’energia, G, in funzione del carico 
applicato, P, per due fissate ampiezze della zona delaminata corrispondenti, rispettivamente, ad a = 0.2 
ed æ = 0.8. A titolo di confronto, le figure riportano sia i valori calcolati con l’espressione (18) sia 
quelli che si ottengono dal più semplice modello a doppia mensola di letteratura, 


=- (21) 


Come si può vedere, i due valori di G tendono ad avvicinarsi al crescere di c, cioè al diminuire 
della lunghezza del tratto dotato del supporto elastico. 
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3(N = , 
Gi) EN | GaN) cel) 
1 - 


a Pon” ; y y o S 
a) b) 
Figura 4 — Velocità di rilascio dell’energia in funzione del carico applicato. 


La figura 5 mostra un confronto tra il carico Pep calcolato con la (20) e quello di letteratura, 


¡Go BEI ey 


Peruerr = pe . 


Anche in questo caso, le previsioni dei due modelli tendono a coincidere al crescere di «. Inoltre, si 
noti come il modello proposto fornisca un valore finito di Pcr anche per & = 0, a differenza di PcrLETT 
che tende ad infinito. In altri termini, è possibile stimare anche il carico di “prima delaminazione” per 
un provino integro, sebbene la validità delle previsioni teoriche così ottenute necessiti di un confronto 
con prove sperimentali. Nell’esempio suddetto, per a = 0, si ottiene Per = 1632 N. 


Figura 5 — Carico critico in funzione dell’ampiezza della zona delaminata. 
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Le figure 6a e 6b riportano i grafici della velocità di rilascio dell’energia, G, in funzione, 
rispettivamente, del carico P ad e; costante e dell’ampiezza di delaminazione adimensionale œ a P 
costante. Dai grafici si può vedere come, al crescere di er, il carico critico Pcr diminuisca mentre, al 
crescere del carico applicato P, il tasso di rilascio dell’energia G aumenti e, di conseguenza, diminuisca 
il valore critico di & per il quale si prevede la propagazione della delaminazione. 


Pe300N 


GM el P= 200 N 4 


a ES ' 03 a 
a) b) 
Figura 6 — Velocita di rilascio dell’energia in funzione del carico e dell’ampiezza di delaminazione. 


A 

250 236.398 | 
a 

200- 

Z 

1504 A - ] 

100 

50- ] 
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wr" e :# 


Figura 7 — Confronto con Laksimi et al. [1990]. 


Il grafico della figura 7 riassume l’intera storia di carico-deformazione di un provino DCB 
sottoposto ad una prova sperimentale a spostamento imposto, 6 = vc. Si distingue un primo tratto 
pressoché lineare della curva, che va dall’origine al punto A e che rappresenta la risposta elastica del 
provino ancora integro (ma provvisto di una delaminazione iniziale di ampiezza ap); nel punto A il 
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carico attinge il valore critico per l’ampiezza di delaminazione assegnata, cioè corrisponde al 
raggiungimento del valore critico della velocità di rilascio dell’energia; il tratto ABC della curva 
rappresenta, quindi, la fase in cui la delaminazione si propaga e il carico diminuisce fino alla completa 
rottura. La figura mostra il confronto tra i risultati ottenuti da Laksimi et al. [1990] (linea grossa) con 


quelli del modello proposto (linea sottile). Si rileva un’ottima corrispondenza tra le due curve. 


60 r ' r r r 
“ s 40 F 
a . 
Z . oS 
3 € 
5 a 
H 199 20 | arbor idi 
e Experiment (Hashemi et al. 1990) ~ (carbomo-cpossidica) 
Anniyneal 
—— FE amiyw 
0 4 4 4 Fi ru "n È 
i p ; a 0 2 4 6 8 
Displacement (mm 5 (mm) 
a) b) 
Figura 8 — Confronto con Zou et al. [2003] per un provino di carbonio-epossidica. 
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Quale ulteriore raffronto con risultati sperimentali, nelle figure 8a e 8b sono confrontate le 
previsioni del modello proposto con i risultati di Zou et al. [2003], relativi ad un provino di carbonio- 


Figura 9 — Confronto con Zou et al. [2003] per un provino di carbonio-PEEK. 
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epossidica caratterizzato dai seguenti parametri: 


l= 100 mm, B=15mm, H=2 mm, ay = 45 mm, 
E, = 130000 N/mm’, E,=3500 N/mm’, Gcr=0.275 N/mm, k = BEJH = 26250 N/mm’. 


Nelle figure 9a e 9b, invece, si riporta un confronto con un secondo caso tratto dal lavoro citato, 
relativo ad un provino di carbonio-PEEK caratterizzato da: 


l= 100 mm, B=15mm, H=2.55 mm, ay = 20 mm, 
E, = 125000 N/mm?, E,=3600 N/mm’, Ger = 2.39 N/mm, k = BEJH = 21180 N/mm’. 


4. CONCLUSIONI 

Nella memoria è stato presentato un modello meccanico del provino DCB, del tipo correntemente 
utilizzato per determinare la resistenza a frattura in modo I dei laminati compositi fibro-rinforzati. Il 
modello proposto consente di prevedere l’inizio e la crescita della delaminazione in condizioni sia di 
carico imposto sia di spostamento controllato. In particolare, grazie all’introduzione di un’interfaccia 
elastica, che tiene conto dell’elasticità della matrice, è possibile calcolare la deformata del provino, 
l’energia potenziale totale e la velocità di rilascio dell’energia. Nel caso limite in cui si trascuri sia 
l’elasticità dell’interfaccia sia quella dei sublaminati, le previsioni del modello coincidono con quelle 
del semplice modello a doppia mensola di letteratura, di cui il presente costituisce una significativa 
estensione. Inoltre, il modello proposto consente di ricavare i valori del carico o dello spostamento 
imposto responsabili dell’innesco e del progredire della delaminazione. La soluzione analitica, ottenuta 
attraverso calcoli lunghi e discretamente complessi, ha ripagato pienamente della fatica spesa per 
ricavarla. Infatti, i primi confronti con dati sperimentali e previsioni teoriche di altri modelli presenti in 
letteratura appaiono molto buoni non solo in termini qualitativi ma anche quantitativi. 
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SOMMARIO: La regione dell’occhio compresa tra il cristallino e la retina contiene il vitreo 
oculare, una sostanza dal comportamento reologico viscoelastico. Nei soggetti giovani il vitreo ha 
la consistenza di un gel ma con gli anni esso perde frequentemente le sue caratteristiche elastiche 
presentando un processo di liquefazione che può estendersi all’intera camera vitreale. Durante i mo- 
vimenti saccadici dell’occhio il vitreo esercita sulla retina tensioni che sono state associate alla pa- 
togenesi del distacco di retina. Nel presente contributo si presentano i primi risultati di un’indagine 
sperimentale sulla dinamica del vitreo oculare realizzata mediante un modello della camera vitreale 
consistente in un cilindro di Perspex recante una cavita riempita con glicerina. Sono stati impiegati 
due diversi modelli: un primo modello con camera sferica, e un secondo con cavita opportunamen- 
te deformata per simulare la presenza del cristallino nella zona anteriore del bulbo. Il modello è 
movimentato da un motore controllato da calcolatore imponendo una legge temporale prestabilita. 
In questa prima fase della ricerca i movimenti saccadici sono stati simulati con rotazioni periodi- 
che. Visualizzazioni del campo di moto del fluido all’interno della cavità sferica, realizzate con un 
tracciante colorato, hanno evidenziato come il moto avvenga essenzialmente sul piani ortogonali 
all’asse di rotazione. Tramite la tecnica PIV sono stati misurati i campi di moto sul piano equatoriale 
ortogonale all’asse di rotazione. Nel caso di cavità sferica i risultati sperimentali risultano in buon 
accordo con risultati analitici basati sul modello semplificato proposto da David et al. [1998]. Si è 
inoltre mostrato come la non sfericità del dominio induca campi di moto dalla struttura assai più 
complessa; in particolare si assiste invariabilmente alla formazione di un vortice, in prossimità della 
lente. 


1. INTRODUZIONE 

L’umor vitreo è la sostanza contenuta all’interno della camera vitreale, la regione dell’occhio 
compresa tra il cristallino e la retina (figura 1a). Dettagliate indagini sperimentali condotte da Lee et 
al. [1992] hanno evidenziato come, dal punto di vista meccanico, il vitreo abbia un comportamento 
viscoelastico. Le proprietà dell’umor vitreo sono, tuttavia, tipicamente molto variabili con l’età del 
soggetto e spesso, si assiste con gli anni ad un progressivo processo di liquefazione, a seguito del 
quale si ha la parziale o totale perdita delle proprietà elastiche del vitreo. 

Esistono nella letteratura medica indicazioni del fatto che la fluodinamica all’interno della ca- 
mera vitreale indotta dai movimenti del bulbo sia legata al processo patologico del distacco di re- 
tina. In particolare, qualora si abbia una rottura della retina, la progressiva infiltrazione del vitreo 
liquefatto attraverso la rottura può portare allo scollamento della retina dall’epitelio pigmentato, lo 
strato ad essa immediatamente più esterno. Tale processo, che viene denominato “distacco di retina 
regmatogeno”, è assai più frequente negli occhi fortemente miopi. 
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Le precedenti osservazioni hanno indotto diversi ricercatori nell’ambito medico a studiare spe- 
rimentalmente tale fenomeno. In particolare è opportuno ricordare il pioneristico lavoro di Lindner 
[1933], che è stato successivamente ripreso da Rosengren e Ostrelin [1976]. Gli autori hanno effet- 
tuato prove sperimentali, su geometria cilindrica, volte a simulare il moto del vitreo a seguito delle 
rotazioni del bulbo oculare. Entrambi i lavori evidenziano l’importanza di ottenere una dettagliata 
conoscenza del moto del vitreo per la comprensione del meccanismo del distacco retinico, seppure 
non permettano alcuna valutazione quantitativa del fenomeno. Il problema è stato successivamente 
riconsiderato da David et al. [1998] che, tramite un’indagine analitica e numerica, hanno superato 
molte delle limitazioni dei precedenti lavori. In particolare gli autori hanno modellato la camera 
vitreale come una sfera posta in rotazione periodica ed il vitreo come una sostanza viscoelastica con 
le caratteristiche reologiche rilevate da Lee et al. [1992]. I risultati degli autori hanno mostrato come 
il moto all’interno del globo risulti sfasato rispetto al moto della parete. Gli autori hanno inoltre evi- 
denziato come le tensioni che il fluido trasmette al contorno aumentino significativamente all’au- 
mentare del raggio della sfera. Tale osservazione è stata interpretata come una possibile spiegazione 
del fatto che gli occhi miopi, tipicamente caratterizzati da dimensioni superiori, sono più frequente- 
mente interessati dal distacco di retina. Gli autori hanno inoltre rilevato come la componente elasti- 
ca del comportamento reologico del vitreo abbia un ruolo modesto rispetto a quella viscosa, sia sul 
moto del fluido che, soprattutto, sulle tensioni che esso trasmette al contorno. 

Recentemente Repetto [2005] ha affrontato il problema del moto del vitreo nella camera vitreale 
mettendo in conto la reale geometria della stessa. Come evidenziato della figura la, infatti, la pre- 
senza, nella parte anteriore del bulbo, del cristallino è responsabile di una rientranza che modifica 
significativamente la geometria sferica. L’autore ha trattato in problema analiticamente nel limite di 
fluido a bassa viscosità. Tale ipotesi risulta appropriata qualora il vitreo sia completamente liquefat- 
to o sia stato sostituito, a seguito di un intervento chirurgico di vitrectomia, con fluidi tamponanti a 
modesta viscosità (spesso a seguito della rimozione del vitreo si pompa all’interno della camera vi- 
treale aria, che viene poi naturalmente sostituita dall’umor acqueo, una sostanza dalle caratteristiche 
meccaniche assimilabili a quelle dell’acqua). L’autore ha messo in evidenza come la non sfericità 
del bulbo induca campi di moto significativamente diversi da quelli che si realizzano in un fluido 
viscoso contenuto in una sfera. In particolare, in prossimità della lente la velocità risulta diretta in 
direzione opposta a quella del moto istantaneo della parete e questo comporta, presumibilmente, la 
generazione di un vortice in tale regione. 


(a) (b) 


Figura 1. (a) schizzo di una sezione verticale dell’occhio umano. 
(b) rappresentazione tridimensionale del dominio utilizzato 
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Nel presente contributo si presentano i primi risultati di un’indagine sperimentale sul moto del 
vitreo all’interno della camera vitreale. Sono stati utilizzati due diversi modelli del bulbo: nel pri- 
mo la geometria considerata è sferica, nel secondo il dominio è stato opportunamente deformato 
per mettere in conto la rientranza dovuta alla presenza delle lente. In figura 1b è riportata una 
rappresentazione tridimensionale del dominio considerato. I movimenti saccadici, ovvero le rapide 
rotazioni del bulbo di ridirezionamento dell’asse visivo, sono stati schematizzati, in via semplificata 
e seguendo l’approccio di David et a/. [1998], come movimenti periodici di opportuna ampiezza e 
frequenza. 

Il presente articolo è organizzato come segue. Nel paragrafo 2 viene decritto l’apparato speri- 
mentale. Il capitolo 3 è dedicato alla descrizione dei risultati conseguiti in geometria sferica mentre 
i risultati relativi al caso di dominio deformato sono riportati nel paragrafo 4. Segue infine un para- 
grafo conclusivo. 


2. DESCRIZIONE DELL’APPARATO SPERIMENTALE 

Il modello di camera vitreale utilizzato consiste in un cilindro in perspex suddiviso in due parti 
uguali, all’interno del quale è ricavata una cavità. Sono state effettuate prove sia utilizzando cavità 
di forma sferica che di forma sferica opportunamente deformata per mettere in conto la reale geo- 
metria della camera vitreale (figura 1b). Il raggio esterno del cilindro è di 120 mm, le dimensioni 
della cavità interna sono di circa 40 mm. La cavità ricavata nel cilindro viene riempita con glicerina 
pura al 98%, ovvero con un fluido Newtoniano di elevata viscosità. Questo significa che, in questa 
prima fase della ricerca, vengono trascurati gli effetti di viscoelasticità del vitreo. Si noti tuttavia che 
il lavoro teorico di David et al. [1998] ha posto in luce come la componente elastica del compor- 
tamento del vitreo abbai un ruolo modesto rispetto a quella viscosa. Essendo l’indice di rifrazione 
della glicerina pressoché identico a quello del Perspex sono stati evitati effetti di distorsione delle 
immagini dovute alla curvatura dell’interfaccia tra le superficie interna del contenitore e il fluido. Il 
cilindro viene montato sul perno di un motore elettrico tipo “brushless” controllato da calcolatore, a 
cui può essere assegnata un’arbitraria legge temporale di rotazione. La posizione istantanea dell’as- 
se del motore viene monitorata tramite controllo analogico. Il motore è sincronizzato con un sistema 
di acquisizione PIV bidimensionale che è stato utilizzato per misurare i campi di moto del fluido. Le 
misure sono state effettuate sul piano diametrale ortogonale all’asse di rotazione. Le particelle trac- 
cianti utilizzate sono spere di vetro cave di diametro approssimativamente pari a 10 um e densità 
pari a quella della glicerina. Uno schema dell’apparato sperimentale è mostrato in figura 2. 


telecamera digitate 


modells de l'occhio 


ta voto ofS 


Figura 2. Schema dell’apparato sperimentale. 
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Ogni campo di velocità (u(x, y, t), v(x, y,£)), con u e v le componenti della velocità rispettivamente 
lungo x e y, è stato misurato 60 volte. Si è quindi potuto valutare un campo (U(x,y, £), V(x,y)) media 
di insieme delle misure. Il corrispondente scarto quadratico medio è sempre risultato inferiore a 
qualche percento dei valori di U e V. La risoluzione spaziale dei campi di velocità ha permesso di 
misurare un vettore velocità per 1.5 mm’. 

La frequenza di acquisizione del sistema PIV è pari a 10 Hz. Al fine di misurare l'andamento 
temporale del moto con maggior dettaglio sono state effettuate acquisizioni multiple, relative a di- 
versi istanti nell’ambito del periodo delle rotazioni; si sono così ottenuti 40 campi di moto per ogni 
periodo. 

L’analisi di cross-correlazione sulle immagini acquisite produce campi bidimensionali in cui i 
vettori sono distribuiti su una maglia cartesiana. Nel caso di esperimenti sul contenitore sferico, i 
dati sono stati successivamente interpolati su una griglia polare (0,4) ottenendo campi di moto in 
termini delle componenti della velocità u, e u,, rispettivamente in direzione radiale e circonferen- 
ziale. Sono stati implementati diversi algoritmi di interpolazione (Agúi e Jemenez, [1987]; Stiier e 
Blaser, [2000]) e a seguito di numerosi test di accuratezza è stato adottato un modello di interpola- 
zione del prim’ordine. Al fine di valutare l’errore sulla stima del singolo vettore velocità interpolato, 
è stato adottata la tecnica Jackknife, in analogia con quanto proposto da Miller [1976]. L’errore sti- 
mato è risultato inferiore a quello commesso sulle velocità medie U e V, garantendo dunque che la 
procedure di interpolazione non ha introdotto errori aggiuntivi ai dati sperimentali. Nel caso, infine, 
di esperimenti su geometria sferica, a seguito della assialsimmetria del dominio, la distribuzione 
radiale della velocità è stata mediata in direzione circonferenziale. Nelle figure 3a,b sono mostrati 
due campi di moto, ottenuti rispettivamente nel caso di dominio sferico e sferico deformato. 

I dati relativi a tutti gli esperimenti effettuati sono riportati in tabella 1 dove gli esperimenti 
su geometria sferica sono stati identificati con la sigla “sfr” e quelli su geometria deformata con 
“def”. 


Exp. | A[deg] | f [Hz] | v[10*ms?] W Exp. | A[deg] | f[Hz] | v[10* ms”] W 
sfr-1 10 1.00 3.79 5.25 def-1 5 1.00 6.30 3.99 
sfr-2 20 1.00 6.58 3.99 def-2 10 1.00 7.62 3.75 
sfr-3 5 1.00 6.25 4.10 def-3 10 1.50 7.45 4.95 
sfr-4 10 0.50 6.79 2.17 def-4 10 1.25 6.98 4.52 
sfr-5 10 1.25 7.16 4.27 def-5 10 1.75 7.12 3.58 
sfr-6 30 1.00 6.73 3.94 def-6 10 2.00 7.22 5.12 
sfr-7 40 1.00 6.5 4.01 def-7 10 2.50 7.02 6.02 
sfr-8 10 1.50 6.11 5.07 def-8 20 2.50 7.13 5.97 
sfr-9 10 1.75 6.26 5.41 

sfr-10 10 0.75 5.58 379 

sfr-11 20 1.50 6.29 5.00 

sfr-12 10 2.00 6.87 5.52 

sfr-13 10 2.50 6.94 6.14 

sfr-14 5 1.00 7.33 3.79 

sfr-15 10 1.25 7.37 4.21 

sfr-16 20 2.50 7.76 5.80 

sfr-17 10 2.00 6.64 5.61 

sfr-18 10 1.50 7.22 4.66 

sfr-19 10 1.75 7.13 5.01 

sfr-20 20 1.50 7.08 4.71 

sfr-21 10 2.50 6.98 6.12 


Tabella 1. Dati sperimentali relativi alle prove su geometria sferica (sfr) e deformata (def). 
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-0.04 40.02 0 0.02 0.04 


Figura 3. Esempi di campi di moto misurati. (a) dominio sferico, (b) dominio deformato. 


3. RISULTATI: IL CASO DI GEOMETRIA SFERICA 

3.1 Visualizzazioni del campo di moto tridimensionale e studio teorico del moto 

Nel problema considerato il moto è indotto dalla rotazione del dominio. La condizione di aderenza 
impone ovunque una velocità puramente circonferenziale su piani ortogonali all’asse di rotazione. È 
quindi lecito aspettarsi che il moto avvenga principalmente su tali piani. Tuttavia la tridimensionalità 
del dominio e la non linearità delle equazioni che governano il sistema possono generare moto anche 
al di fuori di piani ortogonali all’asse di rotazione. Al fine di verificare l'eventuale esistenza di tali moti 
sono state effettuate visualizzazioni all’interno del dominio immettendo, in prossimità della parete e 
appena al di sopra del piano equatoriale, un tracciante colorato. 


Figura 4. Visualizzazioni del campo di moto in falsi colori. Tempi crescenti da sinistra verso destra e 
dall’alto verso il basso. Tra ogni immagine intercorre un tempo pari a 20 periodi (A=20°, f=2.5 Hz). 
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In figura 4 è mostrata una sequenza di immagini (in scala di grigio) che illustra come, sovrap- 
posto al moto primario su piani orizzontali, esista un moto secondario nella forma di due vortici 
toroidali, rispettivamente localizzati nell’emisfero nord ed in quello sud. Nell’emisfero nord le par- 
ticelle fluide si muovono, mediamente, in direzione centripeta quando si trovano in prossimità del 
piano equatoriale, risalgono in prossimità dell’asse di rotazione e discendono quindi in prossimità 
del contorno. Tali moti secondari, la cui esistenza è stata evidenziata numericamente da David et al 
[1998], sono stati invariabilmente osservati, per ogni frequenza ed ampiezza delle oscillazioni. La 
intensità del moto secondario è tuttavia risultata sempre inferiore di 3-4 ordini di grandezza rispetto 
al moto primario circonferenziale. 

La modesta entità delle componenti della velocità permette di trattare il problema secondo un 
approccio analitico semplificato. Supponendo di utilizzare un sistema di coordinate polari sferiche 
(10,9), con r coordinata radiale, O e y coordinate angolari rispettivamente di elevazione e azimutale, 
il vettore velocità u risulta espresso tramite le componenti (u,wu,,u y: Essendo il dominio sferico il 
moto ha caratteristiche di assialsimmetria, ovvero nessuna variabile dipende dalla coordinata 9. 
Sulla base delle visualizzazioni sperimentali è inoltre lecito trascurare le componenti lungo r e 8 
della velocità, scrivendo dunque u=[0,0,u(7,8)]. Tali ipotesi sono alla base del lavoro di David et al. 
che viene nel seguito brevemente descritto per il caso di fluido Newtoniano. L'equazione di Navier 
e Stokes lungo y assume in questo caso la seguente forma lineare. 


Ou, 1 ð| , Ou, 1 0(. , Ou, Us 
=vl> r == sin 13 
Ot r Or Or r”sin0 00 00 r° sin’ 0 a) 
con le condizioni al contorno 
u, = F(t)sind (r=R), (2) 
u, =0 (r=0), (3) 


dove ¢ rappresenta il tempo, v la viscosita cinematica del fluido ed R il raggio della sfera. Nel pre- 
sente caso consideriamo, in analogia a quanto fatto da David et al. (1998), rotazioni periodiche del 
dominio, imponendo dunque F()=Ucos(wt). La linearità della (1) suggerisce di cercare una soluzio- 
ne a variabili separabili nella forma 


u, =e’ 2e(r)sin@ +c.c. 
$ g(r) (4) 


dove c.c. indica il complesso coniugato. Risulta immediato mostrare che le funzione g ammette la 
soluzione 


A 4 
g(r) = ej (lo), == e (Sab) 


con j, funzione di Bessel sferica di prima specie e 
U 


C= ED" (6) 
2 j (KR) 


Nel seguito 1 risultati teorici ottenuti con il modello sopra descritto verranno comparati con le 
misure sperimentali. 
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3.2 Risultati sperimentali 

Come precedentemente descritto si è utilizzato un modello in scala amplificata della camera vi- 
treale. Al fine di utilizzare i risultati conseguiti valutazioni quantitative del moto del vitreo all’inter- 
no dell’occhio è dunque necessario mettere in conto gli effetti di scala derivanti dalle diverse dimen- 
sioni del dominio. Adimensionalizzando, nell’equazione (1), la velocità con U, il tempo con 1/0 e la 
coordinata radiale con R è facile mostrare come l’unico parametro adimensionale che emerge sia il 
numero di Womersley, definito come 


2 
peut 
v 


che rappresenta il rapporto tra le dimensioni del dominio e lo spessore dello strato limite oscillante 
alla parete. Gli esperimenti sono dunque stati condotti conservando tale parametro tra il prototipo 
ed il modello. In particolare si è assunta una viscosità del vitreo (prototipo) pari a v=1.4 10* (Lee et 
al., [1992]) e un raggio della camera vitreale pari a R=0.012 m. In figura 5 sono riportati i confronti 
tra i risultati sperimentali e le previsioni teoriche in termini del modulo della funzione g(r), espressa 
dalla (5a), adimensionalizzata con il suo valore in R. Per aiutare la leggibilità del grafico non sono 
stati riportati le curve relative a tutti gli esperimenti. Ogni curva corrisponde ad un numero di Wo- 
mersley diverso con valori che vanno da 2.4 dell’esperimento 4 a 6.77 dell’esperimento 21. In tutti 
i casi il confronto con le distribuzioni teoriche risulta essere più che soddisfacente, confermando la 
validità delle semplificazioni introdotte nell’analisi teorica descritta precedentemente. 


0 01 02 03 04 05 06 07 08 09 1 
r/r 


Figura 5. Distribuzioni lungo il raggio della funzione adimensionale |g*|=g(r)/\g(R)| 
per diversi esperimenti. I simboli indicano le misure sperimentali e le linee 
le corrispondenti previsioni teoriche 
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4. RISULTATI: IL CASO DI GEOMETRIA SFERICA DEFORMATA 

Nel momento in cui il dominio utilizzato non è perfettamente sferico non è più possibile inter- 
pretare i dati con un modello teorico semplice come quello descritto in precedenza. Nel paragrafo 
introduttivo si erano menzionati i risultati ottenuti da Repetto [2005] nel caso di geometria deforma- 
ta, come nel presente lavoro, ma nell’ipotesi di fluido in moto irrotazionale. Il campo di moto anche 
se con le ipotesi menzionate risultava essere sostanzialmente diverso dal caso sferico, suggerendo 
la possibilità di formazione di una struttura vorticosa legata a fenomeni di distacco dello strato li- 
mite in corrispondenza della deformazione (ingombro dovuto alla presenza del cristallino). I primi 
risultati sperimentali qui riportati confermano che tale supposizione era corretta, mostrando come, 
all’interno di un periodo di un movimento saccadico sinusoidale, si ha la formazione di una struttura 
vorticosa in vicinanza del massimo ingombro durante ogni semiperiodo e che tale struttura migri 
verso il centro del dominio. A partire dai campi di moto bidimensionali ottenuti dalle misure PIV, 
l’esatta posizione del vortice è stata individuata utilizzando un algoritmo di ricerca dei vortici pro- 
posto da Graftieaux et al. [2001]. Prima di scegliere il suddetto algoritmo sono stati implementati 
diversi metodi di ricerca di strutture vorticose, come ad esempio la valutazione del cosiddetto pa- 
rametro di Okubo-Weiss e dell’autovalore del tensore dei gradienti di velocità denominato swirling 
strength. Tuttavia, entrambi i metodi non si sono dimostrati efficienti nelle misure qui descritte, in 
particolare i due metodi non essendo invarianti rispetto a delle rotazioni rigide del sistema di rife- 
rimento, non riuscivano ad isolare la struttura vorticosa di piccola scala rispetto al moto medio del 
fluido anch’esso prevalentemente di rotazione. Al contrario, il metodo proposto da Graftieaux et al. 
[2001] è risultato essere non influenzato dalla rotazione d’insieme dell’intero dominio, riuscendo ad 
identificare correttamente la posizione del vortice. Tale metodo prevede di valutare per ogni punto P 
del dominio discreto il valore di una funzione scalare I°, definita come segue: 


1 
Tl, =— > sin@ 
ay Laine 


dove S è una superficie rettangolare intorno a P contenente N punti in cui si è misurata la velocità e 
6, è Pangolo tra il vettore velocità nell’ennesimo punto e il raggio vettore che unisce tale punto al 
punto P in cui si valuta I. Tale funzione scalare è adimensionale e assume valori compresi tra —1 e 
1, il centro del vortice è individuato dal massimo locale della funzione T, opportunamente soglia- 
ta per eliminare valori inferiori a 0.9. Un esempio della formazione ed evoluzione spaziale della 
struttura vorticosa rilevata nelle misure sperimentali è riportato in figura 6. In tale figura sui campi 
di moto 2D sono state sovrapposte in scala di grigio le curve di livello della funzione I. La regione 
di massimo, contraddistinta da un colore più scuro, è centrata sul centro del vortice. Nella figura 6 
sono riportati otto campi vettoriali che rappresentano otto diversi istanti all’interno di un periodo 
del movimento saccadico sinusoidale equispaziati nel periodo stesso. Il primo campo corrisponde al 
massimo di velocità circonferenziale della parete del dominio (t=0) con un senso di rotazione oraria. 
Si nota la presenza di un vortice già ben formato in corrispondenza della deformazione e un secondo 
vortice in prossimità del centro di rotazione, il quale è il residuo del vortice formatisi in corrispon- 
denza dell’ingombro il semiperiodo precedente. A partire da questo campo di moto, nei successivi, 
corrispondenti agli istanti t= T/8, T/4, 3/8T, T/2, 5/8T, 3/4T e 7/8T, si individua chiaramente la 
formazione del vortice negli istanti t= 3/8T e 7/8T e il suo migrare verso il centro di rotazione nelle 
altre parti del periodo. 
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Figura 6. Campi di moto 2D con sovrapposta la distribuzione delle curve di livello della funzione I al va- 
riare del tempo all’interno di un periodo T dell esperimento def-8. a) t=0; b) t=T/8; c) t=T/4; d) t=3/8T; e) 
t=T/2; f) t=5/8T; g) t=3/4T; h) t=7/8T. 


Gli istanti della formazione corrispondono ai momenti successivi all’inversione di velocità cir- 
conferenziale. In realtà, ci vuole un certo tempo dopo tale inversione per il formarsi di una struttura 
vorticosa. Infatti, a partire dalla posizione del periodo corrispondente all’annullarsi della velocità 
della parete, la distribuzione radiale della velocità presenta un’inversione che a partire dalla parete 
interessa strati di fluido sempre più interni. Il momento di distacco del vortice corrisponde all’istan- 
te in cui tale inversione raggiunge il massimo della deformazione. Da quel momento in poi il vortice 
si sposta quasi orizzontalmente verso il centro del dominio per poi perdere la sua intensità. Analogo 
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comportamento è stato rilevato in tutti gli esperimenti condotti con il dominio deformato, indipen- 
dentemente dall’ampiezza dell’oscillazione e dalla frequenza. La formazione del vortice è quindi 
dovuta esclusivamente alla geometria del dominio. 


5. CONCLUSIONI 

Nel presente lavoro sono state presentate misure bidimensionali del campo di moto del vitreo 
oculare sollecitato da oscillazioni periodiche di tipo sinusoidale. Due configurazioni geometriche 
sono state studiate in dettaglio, una geometria regolare sferica ed una geometria deformata, più 
simile alla forma reale della camera vitreale umana. Nel primo caso i risultati sperimentali, in ter- 
mini di profili radiali della velocità circonferenziale, sono stati confrontati con un modello teorico 
semplificato. Il confronto tra teoria ed esperimenti ha mostrato un buon accordo, confermando la 
validità della teoria formulata da David et al. [1998]. Una seconda serie di esperimenti è stata con- 
dotta con la geometria deformata dalla presenza della lente evidenziando alcune caratteristiche pe- 
culiari della dinamica del vitreo oculare in presenza di un dominio non regolare. Infatti, in tutti gli 
esperimenti svolti si è osservata la formazione di una struttura vorticosa in prossimità del massimo 
ingombro della lente, la quale poi tendeva a migrare verso il centro di rotazione del dominio. Questo 
aspetto rappresenta una prima sostanziale differenza con i campi di moto misurati nel caso sferico, 
in cui non si ha la formazione di nessun moto di piccola scala. I risultati preliminari qui presentati 
riguardano unicamente le strutture del campo di moto. Tuttavia, un’analisi più accurata verrà con- 
dotta per capire come altre grandezze dinamiche vengano modificate dalla presenza di una irrego- 
larità del dominio. Particolare attenzione verrà dedicata alla valutazione delle tensioni tangenziali e 
normali sul contorno del dominio, in quanto tale contorno rappresenta la retina dell’ occhio umano. 
Una quantificazioni degli sforzi indotti dal moto del vitreo oculare sulla retina può avere ricadute in 
campo clinico oculistico nella comprensione dell’insorgenza di alcune patologie legate al distacco 
della retina stessa. 


Il presente lavoro è stato cofinanziato dal MIUR nell’ambito del progetto FIRB2001 (RBAU01Z44F005). Gli autori 
ringraziano l'Ing. Giancarlo Cassini per aver contributo alla messa a punto dell apparato sperimentale e il Prof. Giovanni 
Seminara per le utili discussioni. Un particolare ringraziamento al Dott. Andrea Scupola che ha fornito indicazioni sugli 
aspetti clinici del problema. Federica Di Battista, Claudia Chiesura e Elisa Colangeli hanno contribuito alla presente ricer- 
ca durante le loro tesi di laurea. 
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SUMMARY: Numerical simulations of the flow around a square cylinder between two parallel walls 
(blockage ratio equal to 12.5%) are presented. Reynolds numbers, based on the cylinder side length 
and on the maximum inflow velocity, in the range 100--300 are considered. The numerical solver is 
based on PI finite-element discretization of the viscous terms and on a finite-volume treatment of the 
convective ones on unstructured tetrahedrical grids. Numerical finite-volume fluxes are computed 
through the Roe scheme and a MUSCL-type reconstruction, which leads to second-order accuracy 
and to the introduction of a very small amount of numerical viscosity acting only on the highest re- 
solved frequencies. Preconditioning is used to deal with low Mach numbers, in a formulation which 
maintains time consistency, and time advancing is implicit and second-order accurate. Validation is 
performed by comparing the results of 2D simulations, carried out on progressively refined grids, 
with those of 2D simulations in the literature, characterized by different numerical methods and 
grids. 3D simulations are carried out at Re=200 and 300 and the onset of three-dimensionality in 
the wake after an impulsive start-up is analyzed. 


1. INTRODUCTION 

Bluff-bodies may be used in laminar channel flows for enhancing transport and mixing. Moreover, 
when bluff-bodies with separation points fixed by the geometry are considered, this configuration is 
also interesting for experimental devices (e.g. vortex flowmeters). Recently, the wake of unconfined 
square cylinders in uniform flow at low Reynolds numbers has extensively been studied. However, 
only a few studies have investigated the flow around a square cylinder symmetrically positioned in 
laminar flow between parallel walls (Galletti et al. [2004], Breuer et al. [2000], Li and Humphrey 
[1995], Mukhopadhyay et al. [1992], Davis et al. [1984]). These studies are carried out at low 
Reynolds numbers (Re<1200, based on the square side length and on the maximum velocity of the 
incoming Pouiseuille flow) and they are all two-dimensional. In particular, in Breuer et al. [2000] 
well resolved 2D simulations are carried out at Reynolds numbers up to Re = 300 for a blockage 
ratio equal to 12.5%. In the unconfined case it is shown (Robichaux et al. [1999] and Luo et al. 
[2003]) that, for Reynolds numbers around 160 three dimensional disturbances become unstable in 
the cylinder wake. Although the effect of the confining walls on the critical Reynolds number at 
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which three-dimensional instability of the wake occurs is not clear, it might be possible that the 2D 
approximation used in Breuer et al. [2000] is not valid for all the range of considered Reynolds num- 
bers, as also suggested by the authors themselves. The objective of the present work is to carry out 
3D simulations of the same configuration considered in Breuer et al. [2000] at Reynolds numbers for 
which 3D structures might be present in the wake, and to compare our results with the corresponding 
ones obtained under the hypothesis of two-dimensionality of the flow. Thus, 3D effects on the aero- 
dynamic forces can be characterized for this case and the formation and dynamics of the 3D vortical 
structures in the wake can be analyzed. To this purpose, the code AERO is used, which is a 3D 
solver for the Navier-Stokes equations based on a mixed finite-volume/finite-element formulation 
applicable to unstructured grids (Farhat et al. [1999]), briefly described in Section 2. The test-case 
configuration and the main simulation parameters are described in Section 3. The validation of the 
numerical approach and the analysis of the sensitivity to grid refinement is carried out by comparing 
the results of 2D simulations at different Reynolds numbers with those in Breuer et al. [2000] and in 
Galletti et al. [2004] (Section 4). 3D simulations carried out at Re=200 and 300 are then presented 
in Sec. 5; in particular, the main features of the wake dynamics and structure are illustrated. 


2. NUMERICAL INGREDIENTS 

The Navier-Stokes equations for compressible flows have been discretized in space using a mixed 
finite-volume/finite-element method applied to unstructured tetrahedrizations. The adopted scheme 
is vertex centered, i.e. all degrees of freedom are located at the vertices. P1 Galerkin finite elements 
are used to discretize the diffusive terms. 

A dual finite-volume grid is obtained by building a cell C’; around each vertex è through the rule 
of medians. The convective fluxes are discretized on this tessellation, 1.e. in terms of fluxes through 
the common boundaries shared by neighboring cells. 

The Roe scheme (Roe [1981]) represents the basic upwind component for the numerical evalua- 
tion of the convective fluxes 7: 


F (Wi, 4+ FP (Wy, i) _ 


- W; 
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OF (Wi Wy, n) = yP PRIA a (1) 
in which &* (W;, W;, 11) is the numerical approximation of the flux between the ¿-th and the j-th 
cells, W; is the solution vector at the ¿-th node, ri is the outward normal to the cell boundary and 
R(Wi, Wj, n) is the Roe matrix. The matrix P(W;, W;) is the Turkel-type preconditioning term, 
introduced to avoid accuracy problems at low Mach numbers (Guillard and Viozat [1999]). Note 
that, since it only appears in the upwind part of the numerical fluxes, the scheme remains consistent 
in time and can thus be used for unsteady flow simulations. Finally, the parameter “y, which mul- 
tiplies the upwind part of the scheme, permits a direct control of the numerical viscosity, leading 
to a full upwind scheme for +, = 1 and to a centered scheme when y, = 0. The spatial accuracy 
of this scheme is only first order. The MUSCL linear reconstruction method (“Monotone Upwind 
Schemes for Conservation Laws”), introduced by Van Leer [1977], is employed to increase the or- 
der of accuracy of the Roe scheme. This is obtained by expressing the Roe flux as a function of 
the reconstructed values of W at the cell interface: P” (W;;, Wj, Fi), where W is extrapolated 
from the values of W at nodes z and 7. A reconstruction using a combination of different fami- 
lies of approximate gradients (P1-elementwise gradients and nodal gradients evaluated on different 
tetrahedra) is adopted, which allows a numerical dissipation made of sixth-order space derivatives 
to be obtained. The MUSCL reconstruction is described in detail in Camarri et al. [2004], in which 
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Figure 1: Computational domain 


the capabilities of this scheme in concentrating the numerical viscosity effect on a narrow-band of 
high-frequency fluctuations is also discussed. 

An implicit time marching algorithm is used in the present study, based on a second-order time- 
accurate backward difference scheme. A first-order semi-discretization of the jacobians is used 
together with a defect-correction procedure; the resulting scheme is second-order accurate in time. 

More details on the numerical ingredients used in the present work can be found in Farhat et al. 
[1999] and in Camarri et al. [2004]. 


3. TEST-CASE DESCRIPTION AND SIMULATION PARAMETERS 

The flow around a square cylinder symmetrically positioned between two parallel walls is considered 
here; this configuration and the adopted frame of reference are sketched in Fig. 1. The ratio between 
the cylinder edge length £ and the distance between the walls H is L/H = 1/8. The incoming 
flow is a laminar Pouiseuille flow directed in the æ direction and the considered Reynolds numbers, 
based on the maximum velocity of the incoming flow and on L, range between 100 and 300. Steger- 
Warming decomposition is used to impose non-reflective inflow and outflow boundary conditions. 
At the inflow the Pouiseuille flow is assumed to be undisturbed. Periodic boundary conditions are 
imposed in the spanwise direction and no-slip conditions are forced on the cylinder and on the 
parallel walls. Two different computational domains have been used, for carrying out 2D and 3D 
simulations, which differ only for the spanwise extent of the domain. In both cases, with reference 
to Fig. 1, Li, /£ = 12 and Loy/L = 20. For 2D simulations, the spanwise length adopted 
is #/L = 0.6, and it has been systematically checked that the simulated spanwise velocity was 
negligible. For the 3D simulations, the spanwise length of the domain is #/£ = 6. This choice has 
been made following the indications given in Sohankar et al. [1999] and Saha et al. [2003] for the 
numerical study of the three-dimensional wake instabilities of a square cylinder in an open uniform 
flow. No indication is available in the literature for the case here considered. 

Grid convergence tests have been carried out in the 2D simulations using three grids, mainly 
differing for the spatial resolution in proximity of the cylinder. Details of the grids are reported in 
Tab. 1. We remind that the grids are unstructured and made of tetrahedrical elements. 

The grid GR4, used for the 3D simulations, has been built by replicating the grid GR1 (see Tab. 
1) ten times in the spanwise direction. The details of GR4 are also reported in Tab. 1. 

Since we intend to simulate an incompressible flow, the simulations have been carried out by 
assuming that the maximum Mach number of the inflow profile is Af = 0.1. This value allows 
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nodes on the 
grid tiL 
GRI wa 10 250 0.6 
GR2 210 0.6 
GR3 170 0.6 
GR4 250 6.0 


Table 1: Details of the grids used in the simulations 


compressibility effects to be reasonably neglected in the results and do not pose serious problems 
for the numerics. Moreover, the preconditioning described in Sec. 2 is used to increase the accu- 
racy of the results. In our preliminary simulations it was observed that the preconditioner leads to 
a more accurate value of the pressure coefficient near the stagnation points in the upwind cylinder 
face, improving the mean value of the drag coefficient. The time fluctuations of the force coefficients 
were insensitive to the preconditioner. Concerning the numerical viscosity, the upwind parameter 
“y was set to y = 1.0 on the nodes within a distance equal to 0.1L from the cylinder and y = 0.1 
in the rest of the domain. This choice ensured the stability of all the simulations carried out and, at 
the same time, allowed the preconditioner to be particularly effective in the proximity of the cylinder. 


4. VALIDATION 

In order to validate the numerical approach and to perform a grid convergence study, 2D numerical 
simulations have been carried out for Re=100, 180 and 300, on grids GR1, GR2 and GR3. The 
results have been compared with those obtained for the same configuration in Breuer et al. [2000] 
and in Galletti et al. [2004], with different numerical methods and grid resolutions. The main bulk 
parameters characterizing the aerodynamics forces acting on the cylinder are shown in Tab. 2; as can 
be seen the agreement with the results obtained in the literature is satisfactory. Moreover, it may be 
concluded that grid independence has almost been reached, since in all cases the difference between 
the parameter values obtained with the different grids is very low (< 2%, except for A œp, which 
assumes, however, almost negligible values in most cases). Note that for the Strouhal number, for 
which data are available from both the reference works, the scatter between our results is lower than 
that between the data in the literature. 

As an example of the vorticity fields typical of these 2D simulations, Fig. 2 shows the instanta- 
neous vorticity isocontours obtained in the simulation at Re=180 on grid GR1. Itis worth noting the 
vertical position of the vortices is opposite with respect to the one in the classical von Kármán street; 
indeed, clockwise vortices, which form from the upper edge of the cylinder, are located lower in the 
street than the counterclockwise ones. This behavior is typical of bluff-body flows between parallel 
walls, and has been found, for instance, in Galletti et al. [2004], for the same configuration as the 
one considered here, and in Zovatto and Pedrizzetti [2001], for a confined circular cylinder flow. 

Although this peculiar structure of the vortex street, at Re = 100 and 180 the behavior of the 
aerodynamic forces is the one typically found in 2D simulations of bluff body flows; indeed, as 
shown in Fig. 3 for Re = 180, the time variation of the lift coefficient is perfectly periodic at the 
shedding frequency, while the C’p also exhibits a perfectly periodic oscillation at the harmonic fre- 
quency (not shown here for sake of brevity). At Re = 300 the lift coefficient still shows periodic 
oscillations in time but with an amplitude modulation having a period of approximately three shed- 
ding cycles. A similar behavior, which is absent in 3D simulations (see Sec. 5), has been observed 
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Re St Cp Acro Ac, 
GRI 300 | 0.1250 | 1.6205 | 0.4889 | 3.4019 
GR2 300 | 0.1234 | 1.6359 | 0.5151 | 3.4635 
GR3 300 | 0.1237 | 1.6509 | 0.5191 | 3.4807 


Breuer et al. | 300 | 0.1271 | 1.8603 | 0.5081 | 3.3534 

Galletti et al. | 300 | 0.1320 
GRI 180 | 0.1404 | 1.3659 | 0.0463 | 0.9303 

GR3 180 | 0.1388 | 1.3803 | 0.0459 | 0.9423 

Breuer et al. | 180 | 0.1440 | 1.3250 | 0.0490 | 0.9090 

Galletti et al. | 180 | 0.1370 
GRI 100 | 0.1368 | 1.3758 | 0.0065 | 0.3715 

GR3 100 | 0.1362 | 1.3820 | 0.0066 | 0.3638 

Breuer et al. | 100 | 0.1391 | 1.3500 | 0.0077 | 0.3835 

Galletti et al. | 100 | 0.1386 


GR4(3D) 300 | 0.1345 | 1.4596 | 0.0876 | 1.1889 


Table 2: Main bulk coefficients characterizing the aerodynamics forces acting on the cylinder. St 
is the shedding frequency adimensionalized with £ and the maximum velocity of the incoming 
flow, Cp is the time averaged drag coefficient, Ac, and Aç, are the maximum amplitude of the 
oscillations of the drag and lift coefficients respectively. 
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Figure 2: Instantaneous vorticity isocontours obtained in the 2D simulation at Re=180 on grid GR1. 


in other 2D simulations at the same Reynolds number of the flow around unconfined (Sohankar et al. 
[1999]) and confined (Bruneau, private communication) square cylinders, and might be interpreted 
as a signal of the unsuitability of 2D simulations for this Reynolds number. 


5. 3D SIMULATIONS 
3D simulations have been carried on grid GR4 at Re=200 and 300. For the flow around unconfined 
square cylinders, at these Reynolds numbers, 3D structures are present in the wake showing the main 
characteristics of the so called B transition mode (see Luo et al. [2003] and Robichaux et al. [1999]). 
Our simulations are initiated from an impulsive start-up and the transient phase has been recorded 
and analyzed in order to investigate the mechanisms of formation of 3D structures. 

In Figs. 4 and 5 the time behavior of the lift coefficient is shown for Re = 200 and Re = 300 
respectively. The values of the maximum and minimum spanwise velocity in the field are also 
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Figure 3: Time variation of the lift coefficient in 2D simulations at Re=180 and Re=300 on grid GR1 


reported as an indicator of the occurrence of 3D phenomena in the flow. As can be seen, at Re = 
200, first an almost 2D vortex shedding takes place, characterized by large amplitude oscillations 
(60 < £ < 120). 3D phenomena grow more slowly and their effects on the aerodynamic forces 
become visible only for + > 120, with a significant reduction of the oscillations amplitude of € y 
and of the mean value of the drag coefficient (see Tab. 2), due to the loss of coherence of the vortex 
shedding in the spanwise direction. Conversely, at Re = 300, as expected, the 3D instabilities 
grow more rapidly and 3D effects on the aerodynamic forces are already significant when the vortex 
shedding phenomenon begins to take place; thus, one may say that vortex shedding forms as three- 
dimensional. 
Let us analyze now in more details the form of the 3D instabilities and structures. Two different 
instability modes, initially identified in circular cylinder flows, have also been found for uncon- 
fined square cylinders in experiments (Luo et al., [2003]) and in the Floquet instability analysis 
(Robichaux et al., [1999]). The first one, mode A, occurs at lower Re and is characterized by the 
formation of large scale and wavy vortex loops that connect the spanwise von Kármán vortices. On 
the other hand, mode B is characterized by shorter, finer scaled vortex loops. For unconfined square 
cylinders, mode A was found to occur at Re ~ 160 with a spanwise wavelength of 5.2L, while 
mode B at Re = 190 + 200 with a spanwise wavelength of 1.2L. A third instability mode having a 
spanwise wavelength of 2.8L (mode S) was also identified through the Floquet analysis, which was 
not, however, observed in the experiments. In order to visualize the vortex loops typical of 3D insta- 
bility, in Fig. 6 two isosurfaces of the streamwise vorticity, obtained in the simulation at Re = 200, 
are shown at a time, during the transient, at which the spanwise velocity has reached approximately 
the 10% of the maximum inflow velocity. The structures in Fig. 6 are those typical of mode A of 
instability (see, for instance, the visualizations in Robichaux et al., [1999]), which, however, in our 
case have a spanwise wave length of approximately 6D, due to the imposed periodicity at this length 
in the spanwise direction. The predominance of mode A is maintained also at the latest stages of 
the transient available at present in our simulation; thus, it is plausible to conclude that at this Re 
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Figure 4: Time variation of the lift coefficient and of the maximum and minimum of the spanwise 
velocity obtained in the 3D simulation at Re=200. 
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Figure 5: Time variation of the lift coefficient and of the maximum and minimum of the spanwise 
velocity obtained in the 3D simulation at Re=300. 
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Figure 6: Isosurfaces of streamwise vorticity in the wake obtained at Re = 200 and £ = 92.4. The 
red surface corresponds to &, = 0.1, while the blue one to wẹ = —0.1. 


transition to 3D still follows mode A, while for unconfined cylinders the onset of mode B was found 
at Re = 190 + 200. At Re = 300, since the early stages of transition, the situation is more complex, 
as shown in Fig. 7(a), in which two isosurfaces of the streamwise vorticity are reported at a time at 
which the spanwise velocity has reached approximately the 10% of the maximum inflow velocity. 
Indeed, mode-A-type structures tend to break in smaller vortical loops, not showing however yet 
a well defined spanwise length. Later on in the transient (Fig. 7(b)), only these smaller structures 
are practically visible, which have now a much better defined periodical behavior with a spanwise 
wave-length of approximately 1£ and are, thus, clearly related with the instability mode B. These 
structures persist after the end of the transient (see Fig. 8) and in the developed 3D wake connect 
the vortex tubes of the von Kármán street (in black and light gray in Fig. 8). These spanwise vortex 
tubes are in turn corrugated and distorted by the motion induced by the streamwise vortices, as can 
be seen in Fig. 9, showing the isocontours of the modulus of the projection of vorticity on a (y, 2) 
plane in the wake, together with the velocity vectors (© and w components) on the same plane. In- 
deed, the shear layer is clearly deformed by the vortical motion induced by the streamwise vortex 
loops and this is a typical scenario of 3D wakes at moderate Reynolds numbers, both for circular 
and square cylinders (see, for instance, Luo et al. [2003]). 
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(a) 

(b) 
Figure 7: Isosurfaces of streamwise vorticity in the wake obtained at Re = 300 and (a) £ — 33; 
(b) ~ 37. In (a) the red surface and blue surfaces correspond to wẹ = 0.05 and wy = —0.05 
respectively. In (b) the red surface identifies w = 0.4 and the blue one wy = —0.4. 
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Figure 8: Isosurfaces of the streamwise and spanwise vorticity components in the wake obtained 
at Re = 300 after the transient. The black and light gray surfaces correspond to w. = —0.4 and 
«t= = 0.4 respectively. The pink and blue surfaces identify w = 0.4 and w = —0.4 respectively. 
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Figure 9: Velocity vectors in a (y, z) plane at a = 4 and modulus of the projection of the vorticity 
vector on the same plane; Re = 300 after the transient. 


6. CONCLUDING REMARKS 

Numerical simulations of the flow around a square cylinder between parallel walls, with a blockage 
ratio equal to 8, have been presented. Reynolds numbers, based on the maximum velocity of the 
incoming flow and on the cylinder side length, ranging between 100 and 300 have been considered, 
in order to study the onset of three-dimensionality in the wake. 

The employed numerical solver of the Navier-Stokes equations for compressible flows is based 
on a mixed finite-volume/finite-element formulation applicable to unstructured grids for space dis- 
cretization and on an implicit time advancing. The resulting numerical method is second-order 
accurate both in time and space. 

The validation of the numerical set-up has been performed by comparing the results obtained in 
2D simulations carried out at various Reynolds numbers and grid resolutions with those obtained in 
2D simulations of the same configuration, carried out with different numerical methods and grids, 
available in the literature. The agreement with the reference data is generally good and it may be 
concluded that grid convergence has practically been reached in our simulation; for instance, for 
the vortex-shedding Strouhal number, for which data are available from both the reference works, 
the scatter between our results on different grids is always lower than that between the data in the 
literature. 

Then, 3D have been carried out at Re = 200 and 300 and the transient after an impulsive start-up 
has been analyzed. At Re = 200 the vortex shedding form as two-dimensional and 3D phenomena 
grow up slowly, following the vorticity pattern typical of mode A of instability, previously observed 
in experiments and Floquet analysis. Conversely, at Re = 300 3D effects are noticeable already 
when vortex shedding is beginning and the vorticity configuration is more complex since the early 
stages of transition. At the end of the transient, however, the wake structure is the one typical of mode 
B of transition with streamwise vortices of small spanwise wave-length connecting and deforming 


33 


AIMETA 2005 


the von KérmMan vortex tubes. 

Besides characterizing the differences in the transition scenario between confined and uncon- 
fined cylinders, the next issue of the present study is to investigate whether the transient and fully- 
developed stages of motif can be described (and eventually controlled) in the context of a low-order 
model based on proper orthogonal decomposition. 
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SUMMARY: The capability of theoretical prediction of a separated flow can be fundamental in 
many aspects of industrial design. In this paper an investigation on the capability of the solution of 
the Reynolds Averaged Navier - Stokes (RANS) equations (RANS Technology) to properly predict 
the main features of a separated flow is presented. 

The validation of a numerical method will be referred to basic flows. Applications will be pre- 
sented mainly devoted to the flight technology. Hovwecer considerations and results could be ex- 
tended, generally speaking, to any other industrial application. 


1. INTRODUCTION 

It is well known that the use of the RANS Technology has a number of critical aspects (e.g., a 
proper simulation of the real flow properties, the independence from the numerical parameters), 
even with accurate, converged solutions. This is true not only for separated flows. 

However, accuracy and capability of predicting satisfactorily, with a proper detail, a separated 
flow are potential characteristics, that should be verified for any method and for any application. In 
short, very expensive calculations may be inaccurate or unacceptable. 

On the other hand, RANS Technology can give an approssimate, often acceptable, solution from 
a technical point of view to almost any kind of fluid dynamics problem occurring in the research and 
industry applications. A number of commercial computer codes are available, with a large variety 
of applications. 

On the contrary, more accurate (or exact) numerical techniques (e.g., Large Eddy Simulation or 
Direct Numerical Simulation) cannot be used for practical applications. 

In this paper a 2D RANS method developed by the Authors and other researchers at the University 
of Naples “Federico II’ is presented and discussed. It is well known the a 2D approach for turbulent 
flows has a limitation due to the intrinsic 3D nature of the turbulence. But in the Author”s opinion and 
experience the fundamental features of a 2D experiment could be properly captured. 

With respect to the existing general purpose industrial codes, at least when applied to aerospace 
fluid dynamic problems, this method allows a detailed and accurate numerical analysis of a sepa- 
rated flow. 


2. SOME CRITICAL ASPECTS IN THE PREDICTION OF A SEPARATED FLOW BY US- 
ING RANS TECHNOLOGY 

Here we want to point out some aspects occurring when the RANS Technology is used for the 
prediction of a flow that can be separated. On these aspects the main attention has been focused in 
the development of the method and in the applications that are proposed in this paper. 
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The unsteadiness of a separated flow can be only partly reproduced. The unsteadiness of the 
averaged flow field can be reproduced with some accuracy. 

It is Important to note that in many situations it is not known whether the flow (and the numeri- 
cal solution) is steady or not. So may be critical the choice between a time-consistent (Unsteady) 
and a time-non-consistent (Pseudo Unsteady) solver. The use of the latter is less expensive, but is 
inappropriate if the solution is unsteady. 

The role of the transition may be crucial for some important industrial design applications. So 
its modelling that Is needed in a RANS approach is a crucial point. From the point of view of the 
Authors, it seems clear that fixing the transitions point is equivalent to fix a solution that is not 
‘natural’ when the experiments have been performed with a free transition (i.e., the flow around a 
laminar wing). 

Because of that, the guidelines for the development of the method have been the following: 

- a Time-Consistent approach has been used; 

- the turbulence modelling should have the capability to properly predict the transition. 

One of the main properties of our method is the capability to predict the transition regions, 
even in an unsteady flow. These features may be fundamental for a proper prediction of a separated 
flow. 


3. TURBULENCE MODELLING AND NUMERICAL ASPECTS 

FLOSIM is a 2D URANS code developed at the University of Naples. The code has an explicit 
time marching steady solver with central second order numerical scheme. Unsteady solver uses the 
classical dual time stepping method. 

One and two-equations eddy viscosity turbulence models gave been implemented. Two-equa- 
tions models are k-w standard (wilcox) and SST (Menter) and a non linear k-w EASM. One-equa- 
tion models are the Spalart-Allmaras standard (Spalart) and with curvature correction (SARC model) 
due to Spalart and Shur (Shur). The latter seems to have very good performance for high separated 
and low Reynolds flows also when there is a strong unsteadiness . 

The well known question is the capability of the standard turbulence models to make good pre- 
vision of natural transition. 

Turbulence models often tend to produce too much turbulence, so that the prediction of the tran- 
sition may be inaccurate, above all when an extended laminar zone is expected. 

Typically, there are two ways to solve this problem: 

the first one (and more used) is to fix transition in a point on the wall (2D wall case) consequent- 
ly introducing the turbulence downstream only: of course in this case previous indications about the 
transition are needed; 

another way is to try to reduce turbulence production: this is the method used by FLOSIM by 
using Spalart-Shur turbulence model coupled to a very low free stream turbulence boundary condi- 
tion. 

Some turbulence models seem to be not influenced by the asymptotic turbulence level, while 
other models, such as standard Spalart-Allmaras, may have convergence difficulties when the free 
stream boundary condition is changed. 

The Spalart-Shur turbulence model implemented in FLOSIM Is able to give stable and con- 
verged calculation also when a very low turbulence level is set. 


3.1 Transition on a flat plate 

A fundamental test case is relative to classical no pressure gradient flat plate simulation. Here 
the aim is to evaluate the accuracy of the prevision of the transition zone (Caccavale). Skin friction 
coefficient calculated with Spalart-Shur turbulence model is plotted in figure 1, 1t should be noted 
the good agreement with experimental and theoretical data. 
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Figure 1: Skin friction coefficient on a flat plate 
4, RESULTS 


4.1 Bluff bodies 

Flow around bluff bodies are a widely investigated class of separated flows. Numerical simula- 
tions of the unsteady flow around a circular cylinder with and without a stick (one chord length) at- 
tached at the trailing edge. All the simulations have been performed on computational grids having 
the same characteristics, particularly y+ is always around and less than unit. The Spalart-Allmaras 
with curvature correction turbulence model has been used with low asymptotic turbulence, also in 
the laminar Reynolds number regime. This setting is crucial for a right flow prevision in the tran- 
sitional regime too. As indicated by the turbulent viscosity ratio, at Re=2E+5 turbulent boundary 
layer appears, whereas at Re=1e+5 the flows looks laminar. The calculated mean drag coefficients 
are 1.2 and 0.55 respectively; this result matches well with experimental data (white). In figure 2 


Figure 2: Cl time histories at Reynolds numbers 5.E+5 (left) and 1.E+5 
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the lift coefficient time histories are plotted. The latter shows that turbulence makes vortex shed- 
ding much regular compared to that in the transitional regime. 

An interesting analysis is the comparison between the flow around the circular cylinder in clean 
configuration and with a flat plate on the trailing edge. Figure 3 shows the lift coefficient time his- 
tories at Re=2E+5. It's quite evident the damping effect of the flat plate. 
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Figure 3: Comparison of CI time histories with (right) and without stick at the trailing edge 
(Reynolds number 2.E+5) 


Increasing the Reynolds number the damping effect increases as shown in figure 4. 


Figure 4: CI time history with stick at Re=5.E+5 


In figure 5 it is shown the typical decay of the drag coefficient due to the stick: figure is relative 
to simulation at Re=2.E+5. 
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Figure 5: Drag coefficient time histories with and without stick at Re=2.E+5 


In the next figure a contour map of vorticity in the case of cylinder with stick is shown. 


Figure 6: Vorticity contour map (Re=2.E+5) 


4.2 High lift 

Separated flows for aeronautical applications are analyzed here. The first example here dis- 
cussed are concerned to the prevision of stall mechanism. 

The simulation of flow around NACA0012 airfoil at Mach 0.3 and Reynolds number 3E+6 is 
an interesting application as it involves various phenomena in the stall zone (de Nicola). In figure 7 
the lift curve is plotted, the stall angle is 14 deg with a maximum lift coefficient of 1.4. Calculated 
lift curve has a good agreement with experimental data of Harris: the discrepancy is due to the wall 
interference correction (Harris). 
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Figure 7: Lift curve: NACA0012 airfoil at Re=3E+6, Mach=0. 3 


Our attention is now focused on the flow features at incidence angles around the stall. 

Upper surface skin friction coefficients, plotted in figure 8 for various incidence angles close to 
the stall, show that at angles of attack before the stall there is a laminar bubble that tends to disap- 
pear when the separated zone become larger. 


Figure 8: Skin friction coefficient for various angles of attack (upper surface) 


From figure 9, where a Mach number contour map at 12 deg is shown, it can be noted that a small 
supersonic region ending with a shock wave, appears in the leading edge zone. The shock wave also 
tends to disappear at stall incidence (figure 10). A further analysis at a lower Mach number has 
shown that the bubble isn’t yet present so as the shock wave, this may lead to the hypothesis that the 
bubble seems to be shock induced. All these phenomena make the stall mechanism very complex as 
it isn’t so abrupt as for a bubble explosion but it is due to a trailing edge separation. However, the 
bubble and the shock wave may have a strong influence on this stall. 
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Figure 9: Particular of shock wave at high lift incidence 


Figure 10: NACA0012 airfoil: streamlines in post stall conditions 
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Abrupt stall due to a bubble explosion is well predicted in the next case studied. Has been 
performed numerical simulation of NLR two-component airfoil with 2.6% gap and 20 deg flap 
deflected; Reynolds number is 2.51E+6 and Mach number is 0.185. The computational grid counts 
129920 cells structured in 29 blocks (figure 11). 
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Figure 11: Particular of computational grid for NLR 7301 airfoil 


Calculations have been compared with experimental data (Van Den Berg). These data report the 
existence of a laminar bubble in the leading edge region for all the incidences till the stall, with an 
abrupt loss of lift in the poststall. As shown in figure 12 the calculations confirm the presence of 
laminar separation up to the stall and the bubble explosion at stall angle. 


NLR 7301 airbë Re=2.51E16, Mach=0.158 Ao 


Figure 12: Skin friction coefficient for various incidences 


Figure 13, where lift curve is plotted, confirms the good agreement with experimental data also 
concerning the mechanism of abrupt stall. 
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Figure 13: Lift curve NLR 7301 two-component airfoil 


4.3 Low Reynolds number 

Finally a low Reynolds number unsteady flow calculation around NACA0012 airfoil has been 
carried out at 20 deg incidence. Reynolds number is 2000. An interesting feature is that the calcula- 
tion has been performed with turbulence model activated with low free stream turbulence. The flow 
appears everywhere laminar and this confirms the capabilities of the code to work well at very low 
Reynolds number. In figure 14 the contour map of density and vorticity is plotted, the frame display- 
ing the typical vorticity spot leaving up the trailing edge. A similar result has been found by Iollo et 
al. (Iollo). 


Figure 14: NACA 0012 at Reynolds number 2000, a=20° 
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5. CONCLUSION 

In this paper it has been shown that the use of a proper RANS Technology is able to capture the 
main features of a separated flow. Consequently, an explanation of complex phenomena, so as the 
stall of an airfoil, can be given. 

A 2D method for the prediction of flows with natural transition has been set and discussed; its 
application to separated flows has given satisfactory results. 

In the future our main interest will be addressed to the extension to 3D separated flows. 
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Parole chiave: sistemi meccanici a parametri incerti, analisi ad intervalli 


Sommario 

Questo lavoro nasce dal crescente interesse per il trattamento delle incertezze nella meccanica 
delle strutture. Ci si riferisce in particolare alle tecniche possibilistiche basate sulla matematica 
degli intervalli. L’analisi della letteratura corrente mostra che sono tuttora aperte alcune questioni 
di base legate soprattutto al campo di applicazione e alle prestazioni delle tecniche possibilistiche a 
confronto con tecniche più convenzionali di tipo probabilistico ed all’interpretazione fisica del 
campo di soluzioni ottenibili. Nel lavoro si parte da queste considerazioni per presentare il 
modello di incertezza ad intervalli ed il suo utilizzo nel caso dell’analisi diretta ed inversa di 
sistemi meccanici a parametri incerti. 


1 INTRODUZIONE 


Nella meccanica strutturale il concetto di incertezza viene utilizzato in vari contesti: nella 
modellazione, nell’analisi, nella sperimentazione, nella definizione di grandezze caratteristiche, 
nella definizione di livelli di sicurezza e affidabilità. Quando si affronta un problema 
deterministico è necessario scegliere a priori un modello meccanico rappresentativo del problema 
fisico, nella soluzione invece di un problema a parametri incerti bisogna in aggiunta scegliere un 
modello di rappresentazione delle grandezze incerte. Convenzionalmente la seconda scelta cade su 
modelli probabilistici [Ditlevsen ‘81] anche se esistono modelli diversi per la rappresentazione 
delle incertezze. Ad esempio nel caso della logica fuzzy [Zadeh ‘78] i parametri incerti sono 
rappresentati, in maniera equivalente alla teoria probabilistica, mediante distribuzioni dette di 
appartenenza. Sia le distribuzioni probabilistiche che quelle fuzzy necessitano al loro interno di 
punti di riferimento a cui vengono associati rispettivamente valori di massima probabilità o di 
massimo grado di appartenenza. E” di interesse però prendere in considerazione anche altri tipi di 
grandezze incerte concettualmente distinte dalle precedenti perché non associate a distribuzioni, 
ma definite su intervalli all’interno dei quali la grandezza in questione può assumere tutti i 
possibili valori senza che ve ne sia qualcuno privilegiato. Si dirà allora che la classe di 
appartenenza della grandezza è un intervallo di possibilità. Quest’ultimo tipo di grandezze incerte 
vengono nel seguito meglio formalizzate ed utilizzate. 
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Il modello di incertezza più idoneo al trattamento di un problema meccanico viene scelto in 
funzione della natura delle grandezze incerte che caratterizzano il problema. Si distinguono due 
differenti problemi. Nel problema diretto dato il modello meccanico e data l’incertezza nei 
parametri si vuole valutare l’incertezza nella risposta del modello. Nel problema inverso invece 
note le misure incerte e scelto il modello meccanico si vogliono ricercare quei parametri che 
permettono al modello di riprodurre al meglio le grandezze misurate. 

Nel primo caso l’incertezza associabile ai parametri del modello è frutto della scarsa 
conoscenza dei materiali, della geometria, dalle condizioni di vincolo, ecc. è cioè una misura 
dell’ignoranza che si ha nella fase di modellazione. Nel secondo caso invece, soprattutto nei 
problemi inversi che fanno riferimento all’identificazione costitutiva ed all'identificazione 
dinamica [Friswell e Mottershead ‘95], la scarsa ripetibilità delle misure sperimentali fa nascere 
un'effettiva incertezza nella risposta della struttura. 

In definitiva, nella scelta del modello di incertezza si deve tenere presente che per le incertezze 
di modellazione andranno sempre supposte a priori distribuzioni probabilistiche o fuzzy oppure 
intervalli di possibilità. Nel caso invece delle incertezze sperimentali i diversi modelli di incertezza 
saranno più o meno efficaci in funzione della variabilità e del numero di ripetizioni delle misure a 
disposizione. Infatti la caratterizzazione probabilistica è strettamente dipendente dal numero di 
misure, mentre per la caratterizzazione possibilistica ad intervalli è sufficiente conoscere gli 
estremi di variazione delle misure. 

Obiettivo del presente lavoro è di discutere l’utilizzo delle tecniche possibilistiche nell’analisi 
diretta e inversa di strutture a parametri incerti a confronto con la più convenzionale analisi 
probabilistica. Le metodologie qui presentate sono sviluppate nell’ambito della così detta “analisi 
ad intervalli”, la quale prevede per ogni grandezza incerta l’inclusione in intervalli di possibilità. 
Le soluzioni ad intervalli vengono discusse e confrontate con soluzioni probabilistiche in forma 
chiusa oppure ottenute con simulazioni MonteCarlo. Al fine di rendere coerente e generale 
l’esposizione nella sezione due viene introdotta una notazione unificata per i modelli di incertezza 
a confronto. Nella terza sezione vengono introdotti in maniera formalizzata i problemi incerti 
diretto ed inverso. Nella sezione 4 vengono confrontati alcuni risultati relativi al problema diretto 
nel caso statico e nel caso dinamico con riferimento ad un semplice problema meccanico. 


2 MODELLI DI INCERTEZZA 

Viene utilizzata una notazione unificata per le grandezze incerte, qui indicate in lettere 
maiuscole, indipendentemente dal modello di incertezza con cui vengono rappresentate. Un 
generico parametro incerto è quindi indicato con X. Supponendo che X sia definito sull’insieme dei 
numeri reali, l’incertezza di X è definita dall’intervallo in cui è possibile o probabile che X prenda 
valori. 


Modello ad Intervalli 


Per rappresentare X nell’ambito della matematica ad intervalli, si usa la notazione in 
[Moore ‘66] in cui l’intervallo è definito mediante gli estremi inferiore (xing) e superiore (Xsup): 


X =[x%nt:%ap ]: (1 ER xe SXXX] (1) 


In vista dei confronti con modelli di incertezza convenzionali di tipo probabilistico, 
l’intervallo viene riscritto secondo [Hansen *92] in notazione centrale: 
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inf Xsup 


x= “sa = valore centrale; 


La 
X=x,+AX-e, con AX sa raggio; (2) 


e, = [-1,1] intervallo unità. 


L’espressione (1) evidenzia la natura insiemistica della rappresentazione che costituisce la 
base per la definizione delle operazioni tra intervalli [Moore ‘66]. La (2) è invece una espressione 
di comodo in cui l’incertezza associata ad X è misurata dal raggio dell’intervallo. Si osserva che la 
natura dell’incertezza definita da X è possibilistica nel senso che è possibile che ogni x e X 
assuma valori tra Xing € Xsup Senza che nessuno di questi sia favorito rispetto agli altri. 


Modello probabilistico 

Ai fini del confronto e per coerenza con la rappresentazione (2), si adotta un modello 
probabilistico definito dalla distribuzione uniforme per X nel senso che ciascun elemento 
nell’intervallo ha identica probabilità di essere estratto. In questo caso si ha: 


1 


<x<b 

densità: fy =;b-a a 

0 altrove 

i a+b 

media: E[X] = 

b Di 2 
varianza: o} = ee) 

12 


Si noti che assumendo AX=(b-—a)/2, in coerenza al modello possibilistico, si ottiene 
0? = AX? /3 e quindi: 


X=E[X]ty30, 6) 


cosicché l’espressione (3) di X diviene formalmente identica alla (2) quando si ponga 
E[X]=x,e AX = 3o. 


Incertezza relativa 
In generale, per entrambi i modelli presentati, è conveniente esprimere l’incertezza in forma 
adimensionale. Si introduce allora il rapporto: 


Ex = 


aL (4) 
Xo 
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chiamato incertezza relativa in cui AX è la quantità di incertezza mentre xy è un valore di 
riferimento interno a X, che nel caso possibilistico è il valore centrale x,, mentre nel caso 
probabilistico è la media E[X]. 


3 PROBLEMI MECCANICI CON INCERTEZZA 


Per definitezza, in questo lavoro, la discussione del problema posto è svolta in riferimento ad 
un particolare esempio meccanico. Nello specifico, si considera il problema dell’inflessione di una 
trave incastrata soggetta ad un carico concentrato all’estremo libero (Figura 1). Il modello di trave 
è quello di Eulero e si considerano come parametri incerti l’intensità P del carico e la rigidezza 
flessionale B. Si fa notare come questa scelta permetta di evidenziare sia una dipendenza diretta 
(P) che una dipendenza inversa (B) del modello dai parametri incerti. 


| P 


El = B 


7 
Ubi 


VA 


Figura 1 — Mensola con carico P e rigidezza B a valori incerti 


Si nota inoltre che, da un punto di vista ingegneristico, è di interesse considerare 
l’impostazione sia del problema diretto che di quello inverso e per entrambi la soluzione sia nel 
caso statico che in quello dinamico. Si riassumono di seguito le relazioni che governano i problemi 
citati evidenziandone la dipendenza dai parametri incerti scelti. 


Problema diretto 


Dati P e B incerti si vuole calcolare l’incertezza sulla risposta v che definisce lo spostamento 
trasversale della linea d’asse della trave. 


Caso statico: v(z)= Zg (z) (5) 


Caso dinamico: v(o)=P(0)H (0) == (6) 


D(a) 


In ambito lineare, le equazioni (5) e (6) hanno la stessa struttura ed in esse la funzione v 
rappresenta la soluzione al variare rispettivamente della posizione z oppure della frequenza @. Nel 
primo caso la v(z) è ricavata per integrazione diretta dell’equazione differenziale della linea 
elastica flessionale v" = M(z)/B, essendo M(z) = -Pz il momento flettente. Nel caso dell’analisi 
dinamica, la risposta definita dalla (6) nel dominio della frequenza @è calcolata nell’ipotesi che 
l’intensità del carico P vari nel tempo secondo una legge sinusoidale, in tal caso è possibile 
esprimere ancora la soluzione in funzione del rapporto P/B, attraverso la rigidezza dinamica 


D(o) che si ottiene per B=1. 
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Problema inverso 

Allo stesso modo si presenta un’impostazione formalizzata del problema inverso per mostrare 
il ruolo che giocano le incertezze nella formulazione del problema. Si suppone di conoscere le 
ampiezze della funzione v in alcune posizioni, z oppure @', rispettivamente nel caso statico o 
dinamico. 

In questo caso il problema è allora quello di calcolare i valori di P e B, in modo che la risposta 
analitica v calcolata negli stessi punti z oppure œ sia uguale a quella misurata v”. Questo tipo di 
problemi inversi si definiscono di aggiornamento parametrico ed in letteratura [Sorenson ’80, 
Friswell e Mottershead ’95] vengono classicamente formulati nel seguente modo: 


sP pa 
Caso statico: min(v ,—g{z 7 
( ui ) (7) 
Caso dinamico: min de x : n (8) 
B D(a’) 


In termini generali dunque si deve risolvere un problema di minimizzazione del tipo: 
min(v',v(P,5)) (9) 


in cui il funzionale definito dall’operatore di confronto (+) misura la distanza tra v` e v(P,B). 


E” importante notare che nel caso del problema diretto (cfr. equazioni (5) e (6)) i parametri 
incerti P e B si possono considerare dei coefficienti di valore assegnato ed il valore di v ha come 
variabile indipendente rispettivamente z oppure @. Nel caso del problema inverso invece (cfr. 
equazione (9)), sono i parametri incerti P e B che diventano variabili indipendenti del funzionale 
da minimizzare con la soluzione che viene calcolata in ascisse prefissate (2°, œ). Questa 
considerazione diventa importante quando il problema inverso viene riformulato tenendo conto 
delle incertezze [Gabriele *04]. 

Nelle equazioni dalla (5) alla (9) l'incertezza relativa espressa secondo la (4) assume la forma: 


AP 
Es == 
Po 

10 

LAB ui 
B by 


Per i due modelli di incertezza considerati le quantita definite nella (10) sono esplicitate nelle: 


P= p,(lte,)= p, +AP-e, =[p,.P, | 


Intervalli: 
ies (lte,)=b, + AB-e, =[b,,b,] 


11 
P= p,(1£p) = py +AP=E[P]+430, (1D 


Probabilità: 
B=b,(1£s,)=b, +AB=E[B]+v30, 
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Le espressioni (11) evidenziano come, con la notazione adottata, ci sia una corrispondenza 
diretta nella rappresentazione delle grandezze incerte indipendentemente dal modello di incertezza 
adottato. Si osserva inoltre che mentre per gli intervalli i parametri quantitativi (valore centrale e 
raggio) rimangono fissi, per la probabilità i parametri equivalenti (media e deviazione standard) 
variano al variare della distribuzione ipotizzata e del numero di campioni analizzati. 


4 SOLUZIONE DEL PROBLEMA DIRETTO 


4.1 Caso statico 


Si considera l’inflessione v(z) della trave in Figura 1 definita dall’equazione (5), dalla quale è 
evidente che la soluzione incerta dipende dal valore del rapporto P/B. La soluzione ad intervalli 
segue dall’applicazione della regola di divisione [Moore ‘66]: 


P (1+ EE) 
i _b.(-1+5,)) (2) 


(Cut oA le) a2) 


SA ay = Pelente) o(a) 
b (-1+ £3) 


c 


Nell espressione (12) del valore centrale e del raggio di V si evidenzia come la soluzione è 
non lineare nei parametri £p e & . Per basse incertezze nella rigidezza B il termine £% diviene 
trascurabile rispetto all’unità e la soluzione può essere linearizzata quando l'incertezza di P 
domina su quella di B. Lo stesso non accade per basse incertezze nei valori del carico P. 


Per quanto riguarda la soluzione probabilistica si considerano due diverse soluzioni: 


Ost) 
5 P 
ELV |=» =>8(z): (13) 
Hra Ov ov 0 
By _@P|" 2 aly i av) ) l 3 i 
o;=|—||o0.+—||0=—=% (e, +e, ) 
aP)  LoB,) % 18 


La soluzione a) è fornita da una simulazione MonteCarlo eseguita numericamente su un 
campione generato con distribuzione uniforme. Questa soluzione vista la semplicità del caso in 
esame è anti-economica da un punto di vista computazionale, ma rispecchia una comune prassi nel 
calcolo probabilistico e quindi ha una valenza generale. La soluzione in forma chiusa b) è ottenuta 
mediante sviluppo in serie di v arrestato al primo ordine dell’equazione (5) e fornisce l’espressione 


approssimata di media e varianza calcolate sulla funzione approssimante V . Rispetto alla (12) si 
nota che nella soluzione probabilistica approssimata al primo ordine la media stimata E[Y | 


rimane costante, mentre la varianza è non lineare con proporzionalità diretta a (s&p + ep”), questo 
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implica che per piccole incertezze tali che (£p, eg) << 1 la soluzione risulta sottostimata rispetto a 
quella ad intervalli. 


(I) 


(II) > 


(ID) —— 


Figura 2 — Campo di variazione della linea elastica della trave per incertezze relative £p = eg =0.1. 
Confronto tra intervalli e probabilità (I, IL HI = numerosità campione elevata,media, bassa). 


Per ciascuna delle tre soluzioni incerte ad intervalli e probabilistiche le incertezze relative 
della funzione v sono: 


AV + 
Intervalli : Ey =— = ster, 
V, 146,6 
Bo 
MonteCarlo gae o 
a) VO E D] (14) 
Probabilità : 


30, 
b) Approx 1° ordine e! — a Jere? 
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incertezza su v(2) 


Rapporto tra i fattori di incertezza 


A 
y 
y 


oo 


Figura 3 — Superfici di incertezza per variazioni dell’incertezza relativa fino al 0.5 


Si noti che il coefficiente y3 in e; e ej corrisponde all’aver considerato uniforme anche la 


distribuzione risultante su v. Questa assunzione è stata adottata perché corrisponde ad una 
approssimazione generalmente seguita in problemi più complessi di quello in considerazione in cui 
è difficile stimare analiticamente l’effettiva distribuzione probabilistica della soluzione. 

In Figura 2 sono riportati i risultati della soluzione statica per un’incertezza relativa su P e B 
pari a 0.1. I tre diagrammi sono distinti a partire dall’alto in funzione della soluzione MonteCarlo, 
che è calcolata con un numero di campioni rispettivamente equivalente ad elevato, medio e basso. 
La soluzione ad intervalli ottenuta dalla (12), una volta fissata l’incertezza su P e B, rimane 
invariata ed include la soluzione probabilistica per intero. Un’osservazione importante in questo 
caso è che la soluzione ad intervalli (12) è esatta nell’ambito dell’aritmetica degli intervalli 
[Moore ‘66], questo fa sì che non possano esistere soluzioni fisiche al di fuori di questa. 

La soluzione probabilistica varia invece al variare del numero di campioni. Soluzioni ad 
elevato livello di confidenza si ottengono con altrettanto elevato numero di campioni questo fa si 
che in problemi complessi la componente computazionale risulta penalizzante. D’altro canto le 
soluzioni probabilistiche approssimate permettono una riduzione dei tempi computazionali, ma al 
contrario richiedono un onere elevato negli sviluppi analitici. Le soluzioni ad intervalli sono 
computazionalmente assai più vantaggiose ed hanno l’ulteriore vantaggio di includere sempre la 
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soluzione. Un problema di cui non si discute qui è però la dipendenza dall’espressione 
dell’equazione risolutiva che porta in genere a sovrastimare il risultato [Alefeld e Mayer ‘00] che 
può perdere di significato fisico nel caso di elevata incertezza. 

E” interessante notare come al variare dell’incertezza sui parametri varia l’incertezza su v, è 
inoltre interessante mettere a confronto tale variazione rispetto ai diversi modelli di incertezza. Per 
questo in Figura 3 sono raffigurate le superfici dell’incertezza relativa di v nei tre casi considerati 
(14) in funzione incertezze relative £p, eg. Il grafico in alto mostra che le tre superfici sono simili, 


in particolare si nota che la soluzione MonteCarlo (€; ) è più irregolare perché esclusivamente 
E A E 2 peg n i b 
numerica. Essa è comunque ben approssimata dalla soluzione probabilistica al primo ordine (€; ) e 


che entrambe queste si accordano con la soluzione ad intervalli (ep) la quale come detto tende a 
fornire una sovrastima della soluzione, nel caso attuale il limite inferiore della soluzione 
probabilistica. In definitiva quindi la soluzione ad intervalli è paragonabile a quella probabilistica a 
densità uniforme. Il confronto è quantificato nella superficie raffigurata in basso, in cui si 
rappresentano i rapporti di incertezza tra la soluzione ad intervalli e quelle probabilistiche. In 
questo secondo diagramma si vede, per i motivi sopra detti, che la differenza tra le soluzioni 
aumenta all’aumentare dell’incertezza nei parametri P e B. 


Miss) - estremo Ibero Mies) - estremo ibero 


Figura 4 — Soluzione ad intervalli: (sx) equazione (16), (dx) equazione (17) 


4.2 Caso dinamico 


La soluzione dinamica della trave in Figura 1 è data dall'equazione (6) nell’ipotesi di 
comportamento lineare della struttura e di forzante P variabile nel tempo con un andamento di tipo 
armonico. In questo caso è la soluzione nel dominio della frequenza si può esprimere attraverso la 
sovrapposizione modale [Ewins ‘84] , pertanto la funzione di trasferimento H(@) assume la forma: 


(15) 


in cui i vettori Ø, rappresentano le forme modali della trave normalizzate rispetto alla massa e 
li @ sono gli autovalori reali del sistema non smorzato. 
r 


Il comportamento lineare ipotizzato fa si che la dipendenza della soluzione dall’intensità del 
carico P, sia identica al caso statico, si considera quindi l’incertezza solamente nella rigidezza B. 
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La soluzione ad intervalli degli autovalori di un sistema incerto è un problema non banale 
[Shalaby ‘00], per esso vengono impiegati due metodi di calcolo. Il primo ha carattere generale 
[Qiu, Chen et al. ‘95] mentre il secondo è specializzato all’esempio in questione. La soluzione 
probabilistica di confronto viene calcolata in forma numerica con una simulazione MonteCarlo al 
variare di B distribuito uniformemente. Per la sola soluzione ad intervalli si riportano le equazioni 
risolventi il problema degli autovalori incerti. 


Hia) - estremo ibero 


Figura 5 — Soluzione probabilistica da simulazione MonteCarlo 


Il metodo generale proposto [Qiu, Chen et al. ‘95] è formulato per sistemi discreti e fornisce la 
soluzione ad intervalli distinguendo il problema del limite inferiore da quello del limite superiore. 
Tale metodo soffre di una sovrastima eccessiva della soluzione alle basse frequenze; la correzione 
proposta in [Gabriele ‘04] consiste nel considerare l’intervallo di incertezza iniziale come unione 


di sottointervalli ( B = U 2 ) nei quali il problema della sovrastima è mitigato. La soluzione in 


frequenza è quindi data da: 


2=U2 > H(2)=UA'(2) (16) 


Ù 


dove la Q maiuscola indica l’incertezza nel parametro rispetto al valore deterministico æ. Il 
secondo metodo è invece direttamente derivato dal rapporto di Rayleigh specializzato al caso della 
trave in esame nell’ipotesi di invarianza delle forme modali all’interno dell’intervallo di incertezza 
considerato. L’autovalore r-mo può essere quindi calcolato dall’espressione semplificata: 


Z =B(0Ka,) (17) 


essendo K la matrice di rigidezza calcolata per B = 1. 

Si riportano 1 risultati delle funzioni di risposta in frequenza all’estremo libero della trave, 
ottenute per un’incertezza relativa di B pari a 0.1. In Figura 4 sono riportate le soluzioni ad 
intervalli per i due metodi considerati, la figura di sinistra si riferisce alla soluzione ottenuta con 
l’equazione (16), quindi dall’unione di venti sottointervalli, mentre la figura di destra è riferita alla 
soluzione semplificata a singolo passo ottenuta mediante l’equazione (17). 
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In Figura 5 è riportata la soluzione ottenuta dalla simulazione MonteCarlo in cui sono indicati 
sia i risultati delle singole estrazioni, sia la curva media E(/7), sia le curve E(H Jal Bo, assunte 
come limite di incertezza. Si fa notare che il limite inferiore tende a -o in corrispondenza della 
seconda e della terza risonanza, questo accade quando B30, > E[H ] . Lo stesso accade nel caso 
dell’analisi ad intervalli quando la sovrastima è eccessiva rispetto al valore centrale. Nell’esempio 
la soluzione fornita ad intervalli non soffre di questo difetto in quanto la sovrastima nelle 
frequenze è contenuta dall’aver adottato la partizione B =[),2*. Infine la Figura 6 riporta il 
confronto tra le due soluzioni presentate. Le conclusioni sono identiche a quanto discusso nel caso 
statico in particolare si sottolinea che 1 risultati ottenuti nell’ipotesi di probabilità a distribuzione 
uniforme si accordano con 1 risultati dell’analisi ad intervalli e che questa ultima fornisce sempre 


soluzioni inclusive. La capacità inclusiva dell’analisi ad intervalli ne caratterizza la robustezza, in 
quanto garantisce i limiti di incertezza della soluzione. 


Hi es) - estremo libero 


Figura 6 — Confronto delle soluzioni 


5 CONCLUSIONI 


Le considerazioni esposte nell’articolo sono finalizzate all’analisi del trattamento delle 
incertezze in problemi di meccanica strutturale. Sono state esaminate a confronto due impostazioni 
alternative: una, di tipo convenzionale, basata sulla teoria della probabilità, l’altra di recente 
sviluppo basata sull’analisi ad intervalli. La principale differenza tra le due impostazioni consiste 
nella rappresentazione dell’incertezza che nel primo caso avviene tramite distribuzioni che 
privilegiano un parametro posizionale, p. es. la media, mentre nel secondo caso si fa uso di nozioni 
insiemistiche senza privilegiare valori nell’intervallo di definizione del parametro incerto. Ai fini 
del confronto sono stati considerati modelli semplici di incertezza riproposti secondo una 
formulazione unificata conveniente a rendere omogeneo il confronto. In particolare è stata 
introdotta una misura adimensionale dell’incertezza per la valutazione quantitativa della 
dispersione dei risultati. E” stata mostrata l’impostazione formalizzata per problemi di tipo così 
detto diretto ed inverso sia in campo statico che dinamico con dipendenza parametrica a sua volta 
diretta ed inversa. Le valutazioni quantitative sono state presentate in funzione di un semplice 
problema meccanico nel solo caso di problema diretto. I risultati mostrano che in entrambi i casi, 
statico e dinamico, l’impostazione probabilistica e ad intervalli si equivalgono sebbene con alcune 
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importanti precisazioni a vantaggio della seconda. Infatti, a fronte di problemi legati alla 
sovrastima della soluzione, l’analisi ad intervalli ha la proprietà di includere sempre la soluzione 
fisica del problema con onere computazionale irrilevante rispetto al calcolo in probabilità. 
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SUMMARY: The instability characteristics of a turbulent flat plate boundary layer separating 
under a strong adverse pressure gradient are examined. The analysis is based on the data of direct 
numerical simulation. A theoretical model of harmonic perturbations is considered, including the 
contribution of the turbulent part of the flow, to investigate the stability characteristics of the flow. 
The structure of the organized waves is also investigated by means of Proper Orthogonal Decom- 
positions (POD). 


1. INTRODUCTION 

The understanding of the behavior of a flow in a boundary layer and the knowledge of its physi- 
cal implications like the wall shear rate or the heat transfer is in many technical and industrial ap- 
plications often of paramount importance (for example for the design of a turbine blade). These 
physical properties depend heavily on if the flow is laminar, transitional or turbulent and if a separa- 
tion bubble occurs or not. Therefore, Direct Numerical Simulations of a turbulent boundary layer 
separating under a strong adverse pressure gradient have been performed in Herbst (2004) and 
Herbst & Henningson (2005) which provide a basis for the present work. Flow separation limits the 
performance of many technical devices due to severe pressure losses. For this reason flow control 
may be desired to eliminate or reduce the size of the separation bubble. An interesting control ap- 
proach, investigated in Herbst (2004) and Herbst & Henningson (2005), is to use periodic excita- 
tion to eliminate the separation. The turbulent flow over a flat plate with pressure gradient has been 
chosen as a model for separation and its control for the technological applications of interest, e.g. 
in diffusers. 
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Figure 1: The physical and computational domain with random volume force for 
turbulence generation and the periodic volume force for the separation control. 


These previous studies showed that the parameters (frequency and spanwise wavenumber) of 
the most efficient periodic forcing do not correspond to the most unstable mode predicted by the 
local linear stability analysis. This stability analysis was made based on the time-averaged profiles, 
ignoring the contribution of the turbulent part of the flow. One of the objectives of the present work 
is to advance the analysis by including the effects of the random turbulent flow. For this matter, we 
follow the approach of Hussain & Reynolds (1972) who used an eddy viscosity model to calculate 
the Reynolds stress terms caused by presence of organized waves in a turbulent flow. This approach 
was also used by Reau & Tumin (2000) to investigate the characteristics of harmonic perturbations 
in turbulent wakes. 

The present article is based on the master thesis of S. Speer (2004) and S. Deubelbeiss (2005). 
Since new numerical simulations were performed some parts of the analysis were repeated and 
enhanced. 


2. DIRECT NUMERICAL SIMULATIONS 

The flow in a separated boundary layer over a flat plate has been studied by means of numerical 
simulations. Further, a periodic volume forcing has been used to control the separation bubble. The 
computational box (see figure 1) is 450 nondimensional units long including the fringe region with 
a length 50, 50 units high and 24 units wide. A resolution with 480 points in streamwise direction, 
193 points in wall normal direction and 64 points in spanwise direction is used. Here, the reference 
length is the displacement thickness at inlet of the box where the Reynolds number Re=400. For a 
description of the numerical tools the readers are referred to Herbst (2005). 

At the leading edge x=0 a laminar Blasius boundary layer profile is introduced. Downstream at 
position x=10 the laminar-turbulent transition is triggered by a random volume force near the wall 
(see figure 1). When no control is imposed, the separation of the flow occurs at X „7126 


- - — - — -— — 7 
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x 
Figure 2: Freestream velocity distribution at the beginning of the simulation. 


and the reattachment point is at x,,=247 (reattachment length: x =121). The control is performed 
by means of an oscillating volume force F in the wall-normal direction. This volume force decays 
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exponentially from the boundary layer wall and is centered around x; 


Fy = fo exp (2 -Í y) cos(wt) cos(Bz) , (1) 


T— TO 


Tscale 


where f, is the forcing amplitude, wœ the oscillation frequency, B the spanwise wavenumber, x 
parameter controlling the decay of the forcing in x-direction and c a parameter controlling the wall 
normal decay. It has been shown by Herbst & Henningson (2005) that the closer the forcing position 
(for the separation control) is located with respect to the separation point, the more effective it is in 
suppressing the separation. In the calculations presented here, the location of forcing is x,=110, its 
frequency ©=0.09 and its amplitude set to either f=0.1 or 0.01. 

To give a picture of the flow analyzed here, some plots with fundamental values (from the DNS 
calculations in Herbst (2004) and Herbst & Henningson (2005) are shown in figures 3 and 4. 

In figure 3 the mean flow component in streamwise direction averaged in time and spanwise 
direction can be seen. The upper figure shows the unforced case with the separation bubble (white: 
negative streamwise velocity) and the lower one shows the forced case with a forcing frequency 
w=0.09 and a forcing amplitude f=0.1. For the latter one the separation bubble almost disappears. 

In figure 4 the U -velocity is plotted for the unforced and forced cases. It is obvious that the 
range of fluctuations for the forced case is increased and that its maximum occurs around the forc- 
ing position. 


3. STABILITY ANALYSIS 

3.1. Equations 

Here, we present the equations describing the evolution of organized waves in a turbulent flow. 
The derivation follows closely the work of Hussain & Reynolds (1972). 

As we are interested in evolution of coherent structures in a turbulent flow, we introduce the fol- 
lowing decomposition of variables 


ECZ, t) = E(Z) + ECZ, t) + E (2, t) . (2) 


0 50 100 150 200 250 300 350 400 


Figure 3: Mean velocity component in streamwise direction for the unforced and forced case in con- 

tour plots. Contours of mean velocity component in streamwise direction for (a) the unforced and 

(b) the forced case with the forcing frequency w=0.09 and the forcing amplitude f,=0.1 (neg: dark 
grey to pos: white, white dashed line: -0.025, contour spacing 0.05). 
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where t denotes the mean (time-averaged) quantities, “7 the periodic wave and | j corresponds to 


T 
_ 1 2 
E(#, t) Cp É(X,t+7)dr > (3) 
T J-F 
1 7 
T = T. PY = — P P= == < < 
(Elg, t)) E(z,t ) Him = 2 (5, + iTp) , È a 0 p 27 (4) 


the turbulent motion. The time average £ and phase average <E> are defined as 
a È £ 
Here, T, is the cycle duration and q the phase angle. 


Introducing the triple decomposition (2) into the Navier-Stokes equations, phase averaging 
and subtracting the time-averaged equations give the dynamical equations for the organized waves 
which in non-dimensional form are 

Ou; _ OU; _ OU; Op | 1 07%; lo) 


doti Ox; ti Ox; ~ dz; Re 015 dx; Vila) dx; fij» (5) 
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0 50 100 150 200 250 300 350 400 


Figure 4: U -velocity for the (a) unforced and (b) forced case 
(zero: dark grey to pos: white, contour spacing 0.03). 
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with 
fij = (uu) — uju}. (6) 

Here, u, represents the i-th component of the velocity vector, p the pressure, t the time. The 
Reynolds number is defined as Re=U_“L_/v, /, Where U i L_ and v, represent a characteristic 
velocity, length and kinematic viscosity, respectively. Throughout this text, the subscripts 1, 2 and 
3 correspond to the streamwise (x), wall normal (y) and spanwise (z) directions, respectively. The 
velocities in the respective directions will also be denoted by u=u,, v=u, and w=u,. The term ©. 
is the difference between the phase and time averages of the Reynolds stresses of the backgroun 
turbulence and can be regarded as the oscillation of the background Reynolds stresses due to the 
passage of the organized disturbance. Following the work of Hussain & Reynolds (1972), we adopt 
an eddy viscosity model and assume 


Ou; dù; ` ` 1 /ðù; du; 
fig = pat) = 2455, By = 5 2). 7 
rij des >) sic: A” (= +52) (7) 
Here, v, is calculated based on the DNS results defined either as 
ulul Sy k? 
v=ye= TY or ==: 8 
— 2541511 A E (8) 


with k being the turbulent kinetic energy, e the dissipation rate and c =0.09 a model constant. We 
assume the organized waves to be normal modes given as 
E(x, y, z, t) = E(x, ye’, with 0 = cf o(x')da' + Bz—wt) , (9) 
To 

where a and f are stream- and spanwise wavenumbers, respectively, and © the angular frequency 
of the wave. Introducing the ansatz (9) and (7) in equation (5), removing the terms nonlinear in U, 
and assuming weak streamwise dependency of flow variables result to a modified set of Parabolized 
Stability Equations, PSE (see e.g. Bertolotti et al. 1992). 


3.2 Analysis 
Here we present the results of the stability analysis of the unforced separation bubble. To meas- 
ure the size of the disturbance, we choose an integral quantity defined as 


An? = 131? dy. (10) 
0 
with corresponding growth rate as 


0,2 = =Qí + 2 n( Aa) i (11) 
dx 

In figure 5 the growth rates of disturbance with different B are shown. The data are given for 
both v =v‘ and quasi laminar case v=0. It was observed that for x < 140 the eddy viscosity has a 
damping effect on the growth rate of the disturbances. However, for x >140 the effect of eddy vis- 
cosity was found to be the opposite for B > 0.1. The most interesting result is that when the eddy 
viscosity is taken into the account the damping rate decreases with increasing value of B. The op- 
posite is valid for the quasi laminar case. 

In order to make a comparison with the DNS result, we extract the coherent structures corre- 
sponding to the frequency of periodic forcing by a Fourier Transform of the spanwise- and phase- 
averaged data. Figure 6 demonstrates the structure of the organized wave in the separation bubble. 
As can be seen there, the amplitude of the wave decays as it propagates downstream. 
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In figure 7(a), the values of A ? as a function of streamwise position for the low and high ampli- 
tude cases are plotted. For comparison, the linear PSE results are also given. The DNS data show an 
initial growth of the disturbances which is then followed by a long region of decay. The maximum 
amplitude is reached at x=120 which is close to the point of separation in unforced case x=126. As 
can be seen, the PSE results seem to predict the decay of low amplitude case correctly. Surprisingly, 
the quasi laminar results (v =0) fit better to the DNS data. 
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Figure 5: Growth rate based on A? for different values of spanwise 
wavenumber B. (a) v =v 5 (b) v=0, w=0.09. 
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Figure 6: Contour plot of the coherent structure. Contours ranging from -0.225 to 0.225, 
contour spacing 0.05, dashed lines denote neg. values, solid lines denote pos. values. 


In figure 7(b) the phase of the disturbance (high amplitude case) measured at the location of its 
maximum at each streamwise position is shown. Estimating the wavelength of the disturbance as 
a=d®/dx gives a value of 0.22, which corresponds to a phase velocity of c=0.41. This value is close 
to that found in quasi laminar calculations. 

In figure 8 the amplitude of disturbances from DNS and PSE calculations at two different stream- 
wise positions are compared to each other. Although the variation of 4 ? in PSE and low amplitude 
forcing DNS are in close agreement, the shapes of disturbances are different. 
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Figure 8: Amplitude of disturbances as a function of wall-normal coordinate. 
(a) x=102; (b) x=121. 


4. PROPER ORTHOGONAL DECOMPOSITION 

Proper orthogonal decomposition (POD) was introduced by Lumley to identify coherent struc- 
tures of random turbulent flows. The implementation of the POD in the manuscript is based on the 
method of snapshots developed by Sirovich (1987). For a given number M of snapshots of instanta- 
neous velocity fields u,(x,y,z,t ) at discrete times 1,, the eigenvectors and eigenvalues of the discrete 
correlation matrix 
are computed. For a thorough derivation of the equations and a detailed treatment of the subject, see 
e.g. the Holmes et. al. (1996). 


ner ny nz 3 
Aij = 1/M) 5 y 5 UilTk, Yl, Zl, tiJui(Tk, Yl, Zl, t;) (12) 
k=1 l=1 m=1 n=1 
The results for POD applied to the periodically excited turbulent separation bubble are presented 
in this section. POD has been applied to three dimensional velocity fields, but since the flow is ho- 
mogeneous in spanwise direction, every POD mode is associated with a Fourier mode in spanwise 
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direction. Consequently contour plots of a selected x-y-plane are sufficient to show the results of the 
POD. The contour plots belong all to z=0 and are therefore situated in the middle of the domain. The 
dataset which POD was performed on consists of 150 snapshots taken within one forcing-period 
and equidistant in time. 

Figure 9 shows the eigenfunctions associated with the eight highest eigenvalues. The smaller the 
magnitude of the corresponding eigenvalues the less energy do the eigenfunctions contain. Therefore 
we focus on the eigenfunctions corresponding to the highest eigenvalues. The eigenfunction corre- 
sponding to the highest eigenvalue (mode 1) is shown in figure 9(a) and can be associated with the 
baseflow. Mode 2 and 3, shown in figures 9(b) and (c) show almost similar structures, originating 
around (x=110, y=0), the point where the forcing is centred and decaying only slowly while travel- 
ling upwards in the shear layer. Figure 10(a) and (b) shows the eigenvalues. The first eigenvalue is 
four times higher than the second one. Remarkable is that the following eigenvalues occur in pairs, 
the second and third eigenvalues have the same order of magnitude (0.093 and 0.084) and are more 
than a factor two larger than the fourth and the fifth eigenvalue (0.04 and 0.03). Apart from a phase 
shift the projection of the snapshots on the eigenfunctions representing the time dependence of the 
expansion coefficient of the highest modes (figure 10(c) and (d)) is almost the same for the modes 2 
and 3 with the frequency of the original forcing and for the modes 4 and 5 which represent the first 
harmonic. Since the structures displayed by figures 9(b) and (c) and by the figures 9(d) and (e) are al- 
so very similar is it again likely that each pair belongs to one structure that is travelling downstream. 


5. CONCLUSIONS 

A flat plate boundary layer flow under strong adverse pressure gradient causing a separation 
bubble was considered. The data from direct numerical simulations were analysed by means of sta- 
bility theory and Proper Orthogonal Decomposition technique. 
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Figure 9: Contour plots of the eigenfunctions associated with the eight highest eigenvalues (Mode 1- 8) 
in descending order starting in (a). 
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The organized wave caused by periodic forcing was extracted by Fourier transform of the phase- 
averaged data. 

For stability analysis, a theoretical model of harmonic perturbations is considered, including 
the effects of the random turbulent flow. Here, we have used an eddy viscosity approach to model 
the oscillations of the background Reynolds stresses caused by the organized wave. The equations 
were derived using the nonlocal stability theory based on the PSE method. The stability calculations 
seemed to predict the decay rate of the coherent structures correctly, while the shape of disturbances 
were different. It should be mentioned once again that due to large forcing amplitude the nonlinear 
effects are important here and may be able to explain some of the differences. 

The POD technique was also used to extract the different structures of the flow. These results 
show a similarity to the structures found by Fourier transform of phase-averaged data. 


CN S 


.2 


Figure 10: Eigenvalues corresponding to (a) mode 1- 50 (b) mode 2-8. Projection of 
the snapshots of the DNS on (c) modes 2 and 3 and (d) modes 4 and 5. 
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SUMMARY: This work focuses on the applicability of an U-RANS (Unsteady Reynolds Aver- 
aged Navier-Stokes) model to the simulation of flows around bluff bodies. Test configurations, as 
the cylinder and the square in 2D, and the cube in 3D are discussed for validation purposes. The 
flow around a matrix of cubes has been considered as final application. All computations are per- 
formed by the U-ZEN flow solver, a code developed at CIRA. 


1. INTRODUCTION 

The fluid dynamic analyses of bluff bodies with simple geometries, like cylinders, square and 
cube, represent an interesting problem for the CFD community. These applications are often used 
as test cases for the validation of numerical methods. They are interesting not only in the aerospace 
field, where an airfoil or a wing at high incidence behaves as a bluff body, but also for automo- 
tives and other configurations like civil constructions. The numerical solution of flows with a large 
separation region (comparable to the characteristic length of the domain) is a difficult task because 
of the turbulence that influences the whole field. Although the Navier Stokes equations are able to 
resolve the turbulent flows directly (DNS), the time and spatial resolution needed is so high that 
the applications of industrial interest are not practicable. In fact, the direct simulation requires a 
computational cost proportional to Re” (where Re is the Reynolds number of the problem) [1]. 
Nowadays, the most competitive methodology is based on RANS approach. The weakness of this 
method is the closure problem. The unknowns added by the Reynolds averaging, are treated through 
the turbulence models by introducing an inaccuracy degree not well-defined. The Unsteady RANS 
model is defined by resolving the RANS equations through a time accurate integration. The U- 
RANS computations are more expensive with respect to the RANS, but they allow to resolve fluid 
dynamic problems where an unsteady solution is foreseen. In this way, the domain of the solutions 
extends by including the time dependent ones. Nevertheless, the limitations due to the turbulence 
models continue to be relevant. The aim of this work is to show some applications of U-RANS 
models for flows around bluff bodies. Some theoretical considerations on the applicability of this 
approach are discussed. Two-dimensional applications around the circular and square cylinder are 
performed. The flow around a wall mounted cubical obstacle in a channel flow is computed. Finally, 
the flow around a matrix of cubes placed on the floor of a channel is considered. For all numerical 
simulations the U-ZEN code, developed at CIRA, is used. 


2. CONSIDERATIONS ON THE TIME AVERAGING FOR THE U-RANS MODEL 

The presence in the flow of large recirculation regions is characterized by three-dimensional- 
ity, unsteadiness, and sometimes by macroscopic instabilities such as vortex shedding. For such 
flows, the numerical simulation is rather difficult. The turbulence itself is an intrinsically unsteady 
phenomenon. By excluding a direct numerical simulation, where any time and spatial scale is re- 
solved, the methodologies normally employed in an industrial environment are based on average 
operations. Because of the non-linearity ofthe Navier-Stokes equations, the average operation gives 
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new unknowns that makes the system not resolvable. The closure problem is faced by adding new 
equations that constitute the turbulence model. In the RANS approach, the exact solution ul X, t) is 
decomposed into two contributions: 


u(x,t)= u(x )+u'(x,t) (1) 


where the mean solution u(x) is independent by the time. The turbulent fluctuation is represented 
by u' (x, t). The decomposition (1.) is consistent with the following average operator: 


T 
on | (2) 
u(x) = lim 7 J u(x, t)dt 
As a consequence, 
T 
lim = | u(x,t)dt =0 6) 


T=% 


On the basis of this decomposition, the averaged momentum equation yields new unknowns 
contained in the so-called Reynolds stress tensor. By assuming the definition of the average (2.), 
the RANS equations admit only steady solutions. Besides, even if considering the steady problems 
only, the RANS solution is not equivalent to a time average Navier-Stokes solution because of the 
approximate representation of the Reynolds stress tensor: 


ens # u(x) (4) 


By adding the time derivative of the average solution into the RANS equations, the U-RANS 
model is obtained. But, on the basis of the definition (2.), the presence of a time derivative into the 
momentum equation is not consistent. Therefore, in the U-RANS models the average is defined as: 


T 
u x,t)= > fu x,t Jdt (5) 
-T 


The definition (5.) introduces a new parameter, T. There are several theoretical explanations to 
the use of (5.), but no one is satisfactory enough. A way suggested by Wilcox, [2] is based on the 
concept of the time scale separation. The hypothesis is that at least two time scales T, and T, exist, 
such that T, >> T,, and the time scale T must be included between T << T << T,. In this context, the 
U-RANS should be able to simulate flows where slow variations of the mean quantities with a pe- 
riod T, occur, and the turbulence fluctuations are simulated on the time scale T, by the same closure 
models adopted for RANS. Nevertheless, serious objections are raised to this approach, because 
it is well known that the turbulent flows exhibit more than two time scales, and an energetic ex- 
change among the various scales occur. Besides, the models adopted for the RANS closure problem 
are calibrated through experimental data without any definition of the time scale T, yielding some 
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doubts about the applicability of the average (5.) into the equations. Following these considerations, 
it is possible to state that only certain classes of turbulent unsteady flows can be resolved by the 
U-RANS model. Some cases where the U-RANS have been successfully applied are characterized 
by a single dominant frequency in the field, so that the following decomposition of the solution is 
allowed, [3]: 


u(x,t)= u(x)+ a(x,t)+ u'(x,t) (6) 


The term d(x,t) takes into account the coherent variation of the mean solution. The turbulence 
represented by the term u'(x,t)is not resolved, but modeled as in the RANS approach. So, the cases 
that exhibit an explicit forcing term, such as the simulation of control devices, [6], moving bodies, 
or fluid dynamic instabilities like vortex shedding, [4], [5], can be simulated by an U-RANS model. 
The cases where the turbulence unsteadiness does not present any frequency peak need more ap- 
propriate methodologies. 


3. NUMERICAL MODELING 

The numerical simulations presented here are carried out by using the ZEN code, developed at 
CIRA [7]. ZEN is a very robust, flexible and well assessed computational tool for the aerodynamic 
analysis of complex configuration in subsonic, transonic and supersonic flow regimes. It solves 
the compressible RANS equations around complex 3D configurations (including complete aircraft) 
with a multi-block approach. As discrete mathematical model, a cell centered finite-volume method 
is used, with explicitly added artificial dissipation. A multi-grid method is used with the relaxation 
operator consisting of a Runge-Kutta time integration scheme. The solution procedure is based on 
the time marching concept. The time integration is accelerated by local time stepping and implicit 
residual averaging. The time-marching approach obtains the steady-state solution by evolving in a 
pseudo-time from an initial guess until the convergence is achieved. The time accurate solutions are 
achieved by the dual time stepping method. 


4. FLOW AROUND THE CIRCULAR CYLINDER AT RE 3900 

A first example of flow around a bluff body is the cylinder with circular section. It is a well- 
documented case, for that many experimental and numerical data are available.[8] It is known that 
the flow is laminar and steady for Reynolds up to 40. The flow is still laminar for Re<150, but an 
unsteady vortex shedding occurs. By increasing the Reynolds three dimensional effects are found. 
The sub-critical regime is defined between Re 300 and 2x10°. Here, the separation is laminar, and 
the transition to turbulent flow occurs in the wake. The case examined here is Re =3900, based on 
the cylinder diameter d. Although the experiments show three dimensional effects, the simulation 
has been carried out on a two-dimensional domain. The computational grid, shown in Figure 1, has 
a C-topology, with 200x60 cells. The body circumference has been discretized by 72 cells. The size 
of the first cell near the wall is 3x10*d in the normal direction. The flow domain extends from 20d 
upstream tol5d downstream, where far-field conditions are imposed. Two turbulence models have 
been tested, the standard k-e model, and the TNT Kok k-w model [9]. The transition is imposed at 
50% of the circumference. The comparisons are reported in terms of averaged field over one vortex 
shedding period. The time step used is 0.125d/U_. The Strouhal number obtained numerically is 
0.208, whereas the experimental is 0.215. The mean C, is 0.865, in good agreement with the experi- 
mental value (0.98). In Figure 2 the averaged stream lines are shown. In Figure 3, some velocity 
profiles in the wake are shown. In particular, two abscissas, x/d=1.06 e x/d = 3.00, are examined. 
The z-direction is normal to the streamwise direction. The coordinate system origin is located at the 
center of the cylinder. The experimental data (Lourenco[10]), and some numerical results obtained 
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with LES methods (Mahesh [12], Mittal [13] and Kravchenko [14]), are also reported. For the u- 
component, there are no significant discrepancies between the numerical and experimental data 
(Figure 3a and Figure 3c). The U-RANS results, nevertheless, show weaker gradients in the z-direc- 
tion. At x/d=3.00, Figure 3c, the experimental data show higher values for z/d >1, with respect to 
the numerical ones. Some uncertainties on the velocity values for z/d ® 0 occur. Among the two tur- 
bulence models k-e and k-w, no important difference can be seen. The differences are more evident 
for the normal component, as can be seen in Figure 3b and Figure 3d. The numerical-experimental 
discrepancies are clear, see Figure 3b, moreover the U-RANS disagree with LES for z/d < 0.7. 


5. FLOW AROUND A SQUARE CYLINDER AT RE 22000 

Another example of unsteady turbulent flow is represented by the square cylinder test. The Rey- 
nolds number is 22000, where the reference length is the square edge H. The experimental informa- 
tion are reported in [11]. In such conditions, the flow exhibits a turbulent wake with an alternating 
detachment of vortices. The phenomenology is similar to the previous case. But, here, the separa- 
tions are geometrically imposed because of the sharp edges. The computational grid consists of an 
O topology around the body, see Figure 4. The square edge has been discretized with 48 cells. The 
size of the first cell near the wall is 10H. The flow domain extends from 5H upstream to 14H down- 
stream. Far field conditions are imposed at the external boundaries. A fully turbulent simulation has 
been carried out. To initialize the unsteady computation, a preliminary solution has been computed 
with a large time step. Then, an instantaneous asymmetric disturb has been added. After a transitory 
time of about 60 time units, a periodic regime has been found. The time interval is 0.05, that sam- 
ples the vortex shedding period with about 140 steps. The Strouhal number obtained numerically is 
0.142. The experimental value is 0.132. In Figure 5, the averaged stream lines are reported. It can be 
observed that the flow separates immediately on the upper and lower edges. Two counter-rotating 
vortices are present on the rear face. Besides, near the vertices at x/h=0.5, and z/h=+0.5, two sec- 
ondary vortices are formed. In Figure 6, the comparisons of some velocity profiles are shown. The 
major discrepancies occur in the wake center for z/h O 0. Globally, the results are satisfactory. 


6. FLOW AROUND A WALL MOUNTED CUBE AT RE 40000 

The flow around a cubical obstacle has been the objective of many investigations.[4], [16]. The 
experimental works are due to Hussein [17], and Martinuzzi [18]. The channel height is 2H, where 
H is the cube edge. The Reynolds number is 40000, based on the cube edge and the bulk velocity 
U,. The inflow conditions are imposed at 6H upstream. The out flow boundaries are placed 13H 
downstream. The lateral walls are located 5.5H away from the body. A no-slip boundary condition 
is applied on the floor and on the roof. At the lateral walls a slip-solid wall condition has been used. 
The cube edge is discretized with 48 cells. The size of the first cell near the wall is 10“H. Globally, 
the fine level grid contains 2.4 10* cells. ARANS computation has been carried out by using the 
Kok’s TNT k-o% turbulence model. By looking at the stream lines in the center plane (y/H=0), Fig- 
ure 7, it is possible to note a separation region upstream the cube. The reattachment occurs at x/H 
= -0.74, whereas the experimental value is x/H=-1.04. The computed downstream bubble length is 
2.7H while the experimental value is 1.6H. Another separation is present on the top face of the cube. 
This last one is the trace of the three dimensional arch vortex that surrounds the cube. The upstream 
recirculation, instead, is the trace of the horseshoe vortex. In Figure 8, the pressure coefficient on 
the channel floor and roof along the plane y/H=0 is reported. It can be noted that the pressure on 
the floor downstream the cube is higher than that on the roof, as noted by Hussein [17]. This is due 
to the effect of the horseshoe vortex. By the experiments, the maximum pressure difference (C,, 
= Co...) is 1.47, that is underestimated by the computations, AC,=1.31. Some velocity profiles are 
analyzed in Figure 9. At x/H=0.5 the negative peak is underestimated, whereas, the top wall bound- 
ary layer is well captured. At x/H=1.0 the experiments show a reverse flow that is under-predicted 
by the numerical data. At x/H =2.0 the velocity profile in the bubble center is well reproduced. At 
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x/H=4.0 the numerical results show a slower velocity recovery. Globally, the phenomenology is 
reasonable well reproduced. The main discrepancies occur in the reverse flow regions. 

By using the U-RANS model, Iaccarino (2001), captured a vortex shedding starting from the side 
walls. The arch vortex moved downstream in alternating way, increasing the momentum exchange 
in the recirculation region downstream the cube. So, a better prediction of the bubble length with 
respect to the RANS computation can be obtained through the U-RANS model. Here, the U-RANS 
computation has been carried out on the medium grid level. By adding an asymmetrical instanta- 
neous disturb in the wake, the unsteadiness has not been found. The U-RANS solution converges 
toward the RANS one. This is probably due to the insufficient grid resolution. Future investigations 
will examine this aspect on the finest grid level. 


7. FLOW IN A MATRIX OF SURFACE-MOUNTED CUBES AT RE 3800 

One of the test cases proposed at the workshop held at Delf University in 1997 was the flow 
around a matrix of surface mounted cubes [19]. The experimental data are available on the ER- 
COFTAC website. The matrix is made of 25x10 cubes. The cube edge is H=15 mm. The channel 
height is 3.4H. The distances between the rows and columns have a pitch of 4H in both x and y di- 
rections. In the experimental test, one of the cubes was heated at 75°, but in the present simulations 
the thermic effects are not simulated. The bulk velocity is U, = 3.86 m/s. The Reynolds number, 
based on H and U, is 3800. The experimental measurements were taken at the 18" row from the 
inlet, where the flow is fully developed and symmetric. In this condition, the flow is independent by 
the inlet and outlet, and periodic conditions can be used. The numerical computation is carried out 
on a sub-domain of 4H x 4H x 3.4H in x, y and z directions respectively, see Figure 10. The choice 
of the sub-domain can be arbitrary. In this application the axis origin is shown in Figure 11. The 
periodic conditions are assigned at the faces normal to the x direction, and symmetry conditions are 
used for the faces normal to the y direction. The sub-domain mass flow rate is 13.70 x 10° kg s”. 
The computational grid has 10 blocks, with 626688 cells, Figure 13. The size of the first cell near 
the wall is 2x10°H. The cube edge is discretized by using 48 cells. 

By analyzing the solution in the center plane y/H=0, Figure 12, two different regions can be 
observed. The first is characterized by a turbulent boundary layer on the top-wall, not much influ- 
enced by the presence of the obstacles. For 2H < z < 2.5H, the flow reaches the highest velocity 
in streamwise direction. The lower region, for z<1.5H, is characterized by the separations due the 
obstacles, see Figure 13. The velocity profile at x=-0.3H in the plane y/H=0, intersects the separa- 
tion bubble upstream the cube. Downstream the cube, a primary vortex is present. The vortex struc- 
tures are similar to that of the single wall-mounted cube. A steady RANS computation using the 
Kok’s TNT turbulence model has been carried out. The RANS simulation reproduces most of the 
vortical structures discovered by the experimental data. The major discrepancies are evident in the 
recirculation regions inside the bubbles, that are always over-predicted. The overestimation of the 
separation bubble at the leading edge increases the blockage effect, by diverting the flow laterally. 
This explains the overestimation of the velocities in the plane z/H=0.5 at various abscissas, Figure 
14. The RANS computation underpredicts the streamwise component up to y/H =1.5. An U-RANS 
computation has been carried out starting from the RANS solution, but no differences with respect 
to the steady ones can be observed. An asymmetric perturbation has been applied on the two bub- 
bles upstream and downstream the cube, but, after a transient time, unsteadiness has not been ob- 
served in the solution. In [20], a LES investigation was carried out. By analyzing the energy spectra 
for the streamwise, vertical and spanwise velocity components, no peak was discovered. 


8. CONCLUSIONS 

The U-RANS methodology is not able to resolve every kind of unsteadiness. When in the solu- 
tion a well-defined periodicity is present, it is possible to capture the phenomenology. In this paper, 
four test cases have been shown as examples of flows around bluff bodies. The RANS and U-RANS 
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methodologies have been discussed with respect to the specific applications herein. The two-di- 
mensional cases, both with and without sharp edges in the body, show a clear presence of vortex 
shedding, whereas in the 3D cases this does not happen. In the case of a wall mounted single cube, a 
dynamical mechanism of vortex detachment has been discovered by laccarino et al. [4] by using the 
U-RANS model and by other authors using LES [16]. Some further testing is needed to verify our 
results. For the matrix of cubes an analogous phenomenon has not been found neither in literature 
nor in the present work. 
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Figure 3. Circular Cylinder at Reynolds 3900. Averaged velocity profiles in the wake. (a) x/d = 1.06 u 


component. (b) x/d = 1.06 w component. (c) x/d =3.00 u component. (d) x/d = 3.00 w component 
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Figure 4 — Square Cylinder. Computational Grid 
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Figure 6. Square Cylinder. Averaged velocity profiles in the wake. (a) x/d = 0.875 u component. 
(b) x/d = 0.875 w component. (c) x/d =1.875 u component. (d) x/d = 1.875 w component 
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Figure 8 — Wall Mounted Cube. 
Pressure coefficient on the channel floor and 
roof along the plane y/H =0. 


Figure 7 — Wall Mounted Cube. Stream lines 
in the plane y/H =0. 


Figure 9 — Wall Mounted Cube. Velocity profiles. Streamwise component along the plane y/H =0 
at different stations. (computations: solid lines, experiments: symbols) 
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Figure 11 — Side view of the matrix of cubes on the 
plane y 
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Figure 12 — Streamwise velocity p 
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Figure 13 —Matrix of Cubes. Stream lines. 
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Figure 14 — Velocity profiles in the plane z/H=0.5. Streamwise component. 


77 


Analisi dinamica deterministica ed aleatoria di oscillatori che percorrono 
travi su suolo viscoelastico 


Giuseppe Muscolino 
Dipartimento di Ingegneria Civile, Università degli Studi di Messina 
E-mail: muscolin@ingegneria.unime.it 


Alessandro Palmeri 
Dipartimento di Ingegneria Civile, Università degli Studi di Messina 
E-mail: alexpalm@ingegneria.unime.it 


Keywords: Interazione Veicolo-Armamento, Modello di Zener, Oscillatori Viaggianti, Viscoelasticità. 


SOMMARIO: Si presenta una tecnica per l’analisi nel dominio del tempo delle vibrazioni indotte da un 
oscillatore che percorre una trave elastica su suolo viscoelastico. Le equazioni che governano il moto 
dell’oscillatore mobile e della trave modellata al continuo sono derivate indipendentemente ed accoppiate 
imponendo a valle le condizioni di equilibrio e di compatibilità. Contrariamente alle formulazioni usual- 
mente impiegate nella pratica tecnica, che utilizzano valori equivalenti di rigidezza e di dissipazione vi- 
scosa, l’approccio proposto consente di rappresentate il legame costitutivo della fondazione viscoelastica 
attraverso modelli reologici più accurati. La tecnica proposta è applicata al caso di un oscillatore lineare 
che percorre una trave omogenea poggiante su un tappeto di gomma descritto mediante il modello di Ze- 
ner. Nelle analisi si tiene conto anche dell’effetto dell’irregolarità aleatoria nel contatto oscillatore-trave. 


1. INTRODUZIONE 
In virtù della sua rilevanza nell’analisi e nel progetto dell’armamento ferroviario, la risposta dinamica di 
travi su fondazione elastica soggette a carichi viaggianti è stata oggetto di numerosi studi, teorici e speri- 
mentali (Fryba 1996, 1999). Il modello più semplice è quello della forza viaggiante, in cui l’azione del 
veicolo è descritta da un carico concentrato che si muove lungo il binario. Per tener conto dell’interazione 
dinamica tra veicolo ed armamento si può utilizzare il modello dell’oscillatore viaggiante, in cui il veicolo 
è schematizzato come un oscillatore ad 1 GdL. Così facendo si possono descrivere qualitativamente i 
principali fenomeni di interazione veicolo-armamento. Per una più precisa quantificazione, tuttavia, sono 
necessari modelli più sofisticati, in cui i singoli veicoli ferroviari sono schematizzati come strutture a mol- 
ti GdL. 

L’armamento ferroviario, tuttavia, è usualmente modellato in maniera molto semplificata, come una 
trave elastica poggiante su un letto di molle elastiche con smorzamento puramente viscoso. Nella realtà il 
suo comportamento dinamico è molto più complesso, presentando rigidezza e dissipazione dipendenti 
dalla frequenza di vibrazione. E” questo il caso, ad esempio, di alcune tratte della metropolitana di Mila- 
no, in cui un materassino di materiale elastomerico è interposto tra rotaia e piastra di supporto al fine di 
ridurre l’entità delle vibrazioni indotte dal transito dei veicoli. In questo caso l’assunzione di uno smor- 
zamento puramente viscoso risulta inadeguata, ed è richiesto l’uso di modelli reologici più accurati (Bruni 
& Collina, 2000). 
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Nel dominio delle frequenze l’introduzione di modelli viscoelastici non presenta particolari difficoltà, 
e la funzione di risposta in frequenza dell’armamento ferroviario può essere facilmente valutata a partire 
dalla conoscenza delle proprietà meccaniche dell’elastomero. Tuttavia, per studiare i fenomeni di intera- 
zione veicolo-armamento occorre operare nel dominio del tempo, accoppiando le equazioni che governa- 
no l’equilibrio dinamico delle due sottostrutture (veicolo ed armamento) con le equazioni di stato dei 
componenti viscoelastici. 

In questo lavoro si considera il caso più semplice di un oscillatore ad un 1 GdL che si muove su una 
trave elastica su fondazione viscoelastica, descritta dal modello di Zener. Si perviene ad una formulazione 
nello spazio degli stati in cui le variabili di stato sono lo spostamento e la velocità dell’oscillatore viag- 
giante, gli spostamenti e le velocità dei primi modi di vibrare della trave e le variabili interne aggiuntive 
associate al modello viscoelastico della fondazione. Nello studio si tiene anche conto della irregolarità del 
contatto veicolo-binario, descritta come un processo aleatorio gaussiano, stazionario ed a media nulla, di 
assegnato spettro di potenza. 


2. EQUAZIONI DEL MOTO NEL DOMINIO TEMPO-FREQUENZA 
Si consideri il sistema piano mostrato in Figura 1, costituito da un oscillatore ad 1 GdL che percorre una 
trave elastica su fondazione viscoelastica. 

L’oscillatore viaggiante è costituito da una massa m, che si muove lungo la trave con legge oraria 
nota, x= x,(t), ed è collegata ad essa mediante un dispositivo di sospensione lineare descritto dal model- 
lo di Kelvin-Voigt, in cui la molla elastica k, è in parallelo con il dissipatore viscoso c, . 

La trave, omogenea e ad asse rettilineo, è descritta dal modello di Bernoulli-Navier. I parametri che 
intervengono nell’analisi dinamica sono la lunghezza L, , l’area A, , il momento d’inerzia J, , la densità 
di massa p, , il modulo di elasticità Æ, ed il rapporto di smorzamento viscoso ¢, , che si assume uguale 
per tutti i modi di vibrare. 

La fondazione viscoelastica, infine, schematizzata in figura come un letto di molle, è anch’essa omo- 
genea ed è completamente definita dalla rigidezza dinamica, K,(@), che è una funzione complessa della 
frequenza di vibrazione, @ . 


Ko). 


Figura 1. Oscillatore ad 1 GdL che percorre una trave elastica su fondazione viscoelastica. 
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Nel seguito l’oscillatore viaggiante e la trave su fondazione viscoelastica sono visti come due sotto- 
strutture, per le quali è possibile scrivere separatamente le equazioni di equilibrio dinamico. Queste ultime 
sono accoppiate a valle, tenendo conto delle equazioni di equilibrio e di compatibilità nel punto di contatto. 


2.1. Vibrazioni della trave su suolo viscoelastico 
Il moto della trave rappresentata in Figura 2a è governato dalla seguente equazione differenziale alle deri- 
vate parziali: 


Py, 4&4 Uülx,t)+ E, J, u""(x,6)+D,(x,1)= f(x,6)-K;(0)u(x,t) (1) 


in cui il punto e Papice indicano ordinatamente la derivazione rispetto al tempo ¢ ed all’ascissa x, 
f (x,t) e u(x,t) descrivono rispettivamente il campo delle forze esterne ed il campo degli spostamenti, 
D, (x,1), infine, rappresenta la forza di dissipazione viscosa nella trave. L'equazione (1) è scritta conven- 
zionalmente nel cosiddetto dominio misto tempo-frequenza, utile ad esprimere la dipendenza dalla fre- 
quenza di vibrazione del sistema trave-fondazione. 

Proiettando il moto della trave nello spazio modale, il campo di spostamenti si può esprimere come 
combinazione lineare delle coordinate modali: 


u(x,1)=),,4(09,(0 (2) 


dove f(x) è la i -esima forma modale della trave e q,(t) è la corrispondente coordinata modale, retta 
dall’equazione differenziale ordinaria: 


0+25,0,40+0 4 0)= [SAIL 6) 


in cui ¢, è il rapporto di smorzamento viscoso della trave, che si assume uguale per tutti i modi, ed F, (t) 
è la i -esima forza modale associata alla reazione della fondazione viscoelastica, K,(@)u(x,t), epurata 
della parte puramente elastica, K,(0)u(x,t) : 


v(t) 


b) 


1650) 
u(x,t) 
Ag) 
pe Fete i i 50 


7. 


a) 


x 


Figura 2. Trave elastica su fondazione viscoelastica (a); oscillatore viaggiante (b). 
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ELO e) (4) 
pP,A 


b“ b 


E()= 


essendo K,(0) il modulo di equilibrio della fondazione, che ne rappresenta la rigidezza elastica, ossia la 
rigidezza per carichi statici (al limite per @ > 0). 
La forma modale ¢,(x) e la pulsazione naturale @, nelle equazioni (2) e (3) si determinano risolven- 
do il seguente problema agli autovalori: 
Path ga) = 0,40 5 0, = fa + 
A, 


Pb 


(5) 


con la condizione di ortonormalita: 
Ly 
PA |, AAA, 


essendo 6, l'operatore delta di Kronecker. Si osservi che la forma modale ¢,(x) e l’autovalore œ, non 
dipendono dal modulo di equilibrio, mentre la pulsazione naturale @, aumenta con esso. 

L’equazione (3) è stata ottenuta estendendo al caso di strutture continue con dissipazione viscoelasti- 
ca diffusa il metodo recentemente proposto da Palmeri et al. (2004) per l’analisi modale di strutture inte- 
laiate con dispostivi viscoelastici. Tale metodo consente di eseguire l’analisi dinamica modale di strutture 
a molti GdL dotate di dissipazione viscoelastica tenendo conto in maniera coerente, e senza approssima- 
zione alcuna, dell’effettivo modello reologico dei componenti viscoelastici. 

Nella pratica tecnica, al contrario, si adottano metodi di analisi approssimati, come ad esempio il me- 
todo MSE (Modal Strain Energy), originariamente proposto da Johnson & Kienholz (1982). In accordo a 
tale metodo l’effetto della fondazione viscoelastica sull’ i -esimo modo di vibrare è tenuto in conto me- 
diante valori equivalenti della pulsazione naturale e del rapporto di smorzamento. Così facendo, dunque, 
la trave su fondazione viscoelastica è ricondotta ad un sistema classicamente smorzato, in cui gli oscilla- 
tori modali presentano un modello di Kelvin-Voigt equivalente, con dissipazione puramente viscosa. 


2.2. Vibrazioni dell’oscillatore 
L’equilibrio dinamico dell’oscillatore viaggiante in Figura 2b è governato, nella direzione verticale, 
dall’equazione: 


m,V,()+c,V()+k,v() =m, g (6) 


dove i (t) è l’accelerazione assoluta, v(t) e v(t) sono ordinatamente lo spostamento relativo e la velo- 
cità relativa e g è l’accelerazione di gravità. Quando l’oscillatore percorre la trave (0 < x, (t) < L,) 
l’accelerazione assoluta è somma dell’accelerazione relativa, v(t) , dell’accelerazione del punto di contat- 
to con la trave, d° [u(x, (0), t)]/ dt’, e dell’accelerazione dovuta all’irregolarità del contatto oscillatore- 
trave, d° [-r(x, (t))|/ dt°, essendo r(x) la funzione deterministica che descrive l’irregolarità. Quando, 
invece, l’oscillatore si trova all’esterno della trave (x, (t) <0 o x,(t) > £, ) l'accelerazione assoluta coin- 
cide con l’accelerazione relativa. L'equazione (6), dunque, si può riscrivere nella forma: 
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VO)+2¢, o, VO) + @, Vv) =g- (A) 2 fu 09-10) (7) 


dove @,=k,/m, e €, =c, / [2 O, m, | sono rispettivamente la pulsazione naturale ed il rapporto di 
smorzamento dell’ oscillatore viaggiante e dove II, (£) è una funzione finestra che vale 1 quando il vei- 
colo è in contatto con la trave, 0 in caso contrario: 


II, (1) =U (x, (M)U(L, =x, (1) (8) 


essendo U(x) la funzione gradino unitario, tale che U(x)=0 per x<0 e U(x)=1 per x>0. 


2.3. Condizione di compatibilità 
Si osservi che l’equazione (7) garantisce la compatibilità tra le vibrazioni di trave ed oscillatore mobile. 
Esplicitando l’accelerazione relativa a secondo membro, si ottiene: 


W0+2£,0, Wt) +00) = g+ N(t)+ 


1, (0) [(x,(0,0)+24',0,01,(0+4',0,9%(0+4'x,(0,03,(0 | di 


essendo N(#)=II,(t) [Gx (1) n (t)+r'(x, (A) x, e] la forzante dinamica sull’oscillatore viaggiante, 
la quale, nota l’irregolarità, r(x), e la legge oraria, x,(t), può essere valutata preventivamente. Sosti- 
tuendo, quindi, l’equazione (2) nella (9) si ottiene: 


BO +26, 0,00 + OVO = Ft NO- LV {MOG0+C0G0+K, Oa} (10) 


dove i coefficienti tempo-dipendenti M,(t)=¢(x,O)0,@, C(M)=24G (0x0, 0 e KA= 
[da i+ (MET, (0) sono rispettivamente la massa, la dissipazione e la rigidezza con 
cui l’ i -esima forma modale della trave interagisce con le vibrazioni verticali dell’oscillatore viaggiante. 


2.4. Condizione di equilibrio 
Al fine di soddisfare l’equilibrio si esprime la forzante f(x,t) nell’equazione (3) come la reazione del 
dispositivo di sospensione dell’oscillatore viaggiante, localizzata nel punto di contatto con la trave, 
x=x (t): 


S&A [k v(t) wO, OEE- xO) (11) 
dove 6(x) è la funzione delta di Dirac. Sostituendo l’equazione (11) nella (3), si trova: 


G(t)+2¢,0,9,() +o; q,(0=0,O[x,v0+c,40]9(,00)-E,0 (12) 


3. MODELLO DI ZENER 
In Figura 3 é rappresentato il modello di Zener, costituito da una molla k, in serie con un modello solido 
standard, in cui una molla k, € in parallelo con un elemento di Maxwell, ottenuto disponendo una molla 
k, in serie con un dissipatore c,. Tale modello è stato utilmente impiegato da Bruni e Collina (2000) per 
rappresentare il comportamento dinamico del materiale elastomerico che in alcune tratte della metropolita 
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Fo) 


k, 


Figura 3. Modello di Zener. 


di Milano è stato inserito per ridurre gli effetti dinamici su veicoli ed armamento. Si dimostra che la rigi- 
dezza dinamica di questo modello è: 


-1 Al 
K(0)= paz | 
joc, 


ed il modulo di equilibrio è K(0) =(k;'+k,')". 
Nel dominio misto tempo-frequenza la forza di reazione del modello viscoelastico, epurata dalla parte 
puramente elastica, F(t) , € legata alla deformazione, d(t) , dalla relazione: 


-1 


F4) =[K(@)- K(0)]d(1) (13) 
In alternativa, la forza F(t) nel dominio del tempo si può esprimere come: 
FO =k A(0+k, 2, (0 -K(0)d(t) (14) 


essendo A (t) e 4, (t) due variabili interne che misurano la deformazione delle due molle k, e k, . Que- 
ste ultime sono governate dalle equazioni: 


miro LO . _% k 
4,(t)= 40) x A) Li Ran (15) 


dove 7, =k,/c, è il tempo di rilassamento dell’elemento di Maxwell. Sostituendo la seconda delle equa- 
zioni (15) nella (14) e nella prima delle (15) si ottiene: 


FO=a440 5 A4(9=D,d(0)+D,4(0) (16) 
avendo posto: 
k k 
O E yea e 
k, +k, k, +k, +k, (k, +k, +k,)7, 


4. MODELLO ALEATORIO DELL’IRREGOLARITA’ 
L’irregolarità nel contatto oscillatore-trave, r(x), si può pensare come la realizzazione di un processo 
aleatorio, gaussiano, stazionario ed a media nulla. Questo, dunque, è completamente definito dalla spettro 
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di potenza monolatero, G,(y) , dove y è il numero d’onda. Assegnato lo spettro G,(y) , il generico cam- 
pione dell’irregolarità si può generare mediante la formula: 


r(x)=)"_A4;sin(y,x+0,) 


dove A, =,/2G,(y,)Ay e y;=(j-1/2)Ay sono rispettivamente l’ampiezza ed il numero d’onda del 
j -esimo contributo, Ay =y,/n è il passo di campionamento, y, è il numero d’onda di cut-off, n è il 
numero di termini considerati nella sommatoria e 9; è la j -esima fase aleatoria, uniformemente distri- 
buita nell’intervallo [0,27] . Il corrispondente segnale N(t) si può esprimere nella forma: 


M(t) =T, (Ù [-20) Y 47 sing, x+8,) +80) Oo" 4,7, cos(y, x+0,)| (17) 


In ambito ferroviario (Fryba 1996), ad esempio, si può utilizzare per l’irregolarità r(x) il modello 
empirico proposto dalla francese SNCF (Société Nationale des Chemins de Fer): 


_ 10° A, 
(+7/r)° 


dove A, é un parametro delle dimensioni di una lunghezza al cubo, legato alla qualita dell’armamento, e 
y, è una costante, delle dimensioni di un numero d’onda. 


GY) (18) 


5. RISPOSTA NELLO SPAZIO DEGLI STATI 


5.1. Equazioni accoppiate 
Nei precedenti paragrafi si sono presentate le equazioni che governano le coordinate modali della trave 
elastica (eq. (12)), le vibrazioni dell’ oscillatore viaggiante (eq. (10)) ed il modello di Zener (eq. (16)). Tali 
equazioni possono essere riscritte in forma matriciale, una volta introdotti i vettori delle variabili di stato 
z(0=[9,(0 å 0] per l’ i -esimo modo di vibrare della trave e z, (£) =[v(t) vo] per l’oscillatore 
viaggiante. Inoltre, dal confronto tra le equazioni (4), (13) e (16) si trae che una variabile interna, 4, (t), 
i=1,2,---, deve essere aggiunta per ciascun modo di vibrare della trave. 

Così facendo, l’equazione (12) diviene: 


2,(0)=D,2,(0+D,(0z,()-vE() (19) 


in cui: 


v 


p =| ° i ; D= SEB Haas 
E -0 26,0, > (1) =P, (x,(0) v(t) k c, > V= 1 


L’equazione (10) diviene: 
i.0=D,2,(60+v[g+N0]+Y[D,020-vM0v' 2,(0)] (20) 


in cui: 
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l 0 1 | l 0 0 | 
D, = 2 ; D, (t) = 
O, -2 È, O, -K; (t) -C, (t) 


Sostituendo I’ equazione (19) nella (20) si ottiene: 
2,(0=D,(0z,(0+9 D,(0z,(0+v [g+ NO]+v [MOL] (21) 
avendo posto: 


D,()=D,-vv' MOD, ; D,(0=D,()-vv' MOD, 


Tenuto conto della prima delle equazioni (16), le forze F, (t) nelle equazioni (19) e (21) si possono 
esprimere nella forma: 


a 
Ps A, 


F(t) = Ay 0) (22) 


in cui la variabile interna /,,(¢) si può pensare come la deformazione nella molla k, nel modello di Zener 
che descrive l’effetto della fondazione viscoelastica sull’ i -esima forma modale della trave. La variabile 
interna /,,(¢) , in virtù della seconda delle equazioni (16), è governata dall’equazione: 


ÅO =d? Z,(0+D,4,(0) (23) 


dove dj =[0 D,]. 

Le equazioni (19), (21), (22) e (23) devono essere risolte simultaneamente. A tal fine, fissato il nume- 
ro m di modi della trave da considerare nell’analisi, si introduce il super-vettore delle variabili di stato 
H(t)=[ 2) ZO A() += A,(0) ||, che è retto dall’equazione: 


i(t) = DE) z(t) +f (t) (24) 


dove: 


D,,, (1) | E 


[v [2 +N(0)]| 


; f©=| [0 oy 


D(*) = 
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5.2. Schema numerico 
L’equazione (24) rappresenta un sistema di 3m+2 equazioni differenziali lineari ordinarie non omoge- 
neee, di primo grado, con coefficienti tempo-dipendenti, la cui soluzione nel dominio del tempo può esse- 
re conseguita mediante la soluzione al passo descritta qui di seguito. 

Si discretezza l’asse dei tempi in intervalli di ampiezza costante, At. Nell’ i -esimo intervallo tempo- 
rale, [¢,,t,,,], con t, =i At, si assume che la matrice dinamica sia costante, D, = D(*,). Conseguentemen- 
te, all’interno di tale intervallo l’equazione (24) si può pensare nella forma: 


Z(0) =D, z(0)+[D(1)-D,]z()+f(0) , 1,,<t<t, 
dove a secondo membro, oltre al vettore forzante f(t), compare una pseudo-forzante, anch’essa tempo 
dipendente, LO) -D,|z() . Se quest’ultimo termine è lineare nel passo, il vettore della risposta a fine 
passo, Z,,, = Z(t,,,) , si può valutare mediante l’espressione: 


Za =0,Z, +Ya (D,a -D, Jm + Y jo ft) 4+ Ya fa) (25) 
dove le matrici ®,, y,, € Y, sono così definite: 


©, = exp(D, At) ; L, = (0, -Im )D7' > Yio =(0,-1L,)D;' > Ya =( gl -Ln ) D7 


i Ai 


essendo I,,,,, la matrice identità di ordine 3m +2 . Risolvendo l’equazione (25) rispetto a z,,, si trova 


(Muscolino 1996): 


i+] 


Zi = O, Z, + Tio f(t) + Ya fC) (26) 
in cui: 


=1 ~ pe Gai 
J, =[ Ln Ya (D,a -D,)| ; 0,=3,0, 3 Vo=]: Yo 5 Ya=3J,Y5 

Lo schema numerico così ottenuto è stato implementato su software Mathematica 4.0 per eseguire le 
prove numeriche presentate nel paragrafo successivo. Per i confronti, uno schema numerico del tutto si- 
mile è stato formulato ed implementato per il caso in cui la fondazione viscoelastica è descritta dalla me- 
todo MSE. 


6. APPLICAZIONI NUMERICHE 
Al fine di validare la tecnica proposta, e di evidenziare l’errore che si commette utilizzando il metodo 
MSE, si sono studiate le vibrazioni indotte dal passaggio di un oscillatore che percorre una trave sempli- 
cemente appoggiata con velocità costante V , per cui la legge oraria dell’oscillatore viaggiante è 
x,(t)=Vt. Le caratteristiche geometriche e meccaniche di trave, fondazione ed oscillatore viaggiante 
sono riportate in Tabella 1. 

La i -esima forma modale della trave è data dall’ espressione: ¢ (x)= /2/m, sin(i7x/L,), essendo 
m, = Pp A, Ly la massa della trave. In Figura 4 sono rappresentate le prime cinque forme modali, ed in 
Tabella 2 sono riportati i corrispondenti valori delle frequenze circolari e dei rapporti di smorzamento che 
si ottengono utilizzando la tecnica proposta ed il metodo MSE. Si noti che i primi modi di vibrare presen- 
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Tabella 1. Dati geometrici e meccanici. 


trave elastica oscillatore viaggiante fondazione viscoelastica 


Po 8000 kgm | my 120 kg ik 5000. KN/m? 

E, 300000 Mpa i k 10 kNm i ci 100 kNs/m? 

L, 10 m i cy 200 Ns/m i ka 1000 kN/m? 

di 0.02 y 10 ms È k 2000 kN/m? 

Ay 50 cm? o, DS 9.13 rad/s Ma IE è oo 
Jy 417 em È € 0.0913 i k0) 667 kN/m? 


tano uno smorzamento maggiore, legato alla dissipazione di energia aggiuntiva che è garantita dalla fon- 
dazione elastomerica. 

In una prima fase si è trascurato l’effetto della irregolarità del contatto oscillatore-trave, ponendo de- 
terministicamente N(t)=0 nell’equazione del moto. In Figura 5 sono mostrate le deformate della trave 
quando la posizione dell’oscillatore viaggiante, x= x,(t), è 0.25£, (a), 0.50L, (b) e 0.75L, (c). In 
tutte e tre le configurazioni si osserva che il metodo MSE, utilizzando valori equivalenti di rigidezza e 
dissipazione viscosa, sottostima notevolmente le deformazioni subite dalla trave per il passaggio 
dell’oscillatore. 

In Figura 6 sono mostrate le storie temporali dell’abbassamento in mezzeria (a), u, (t) =u(L, /2,t), 
dell’abbassamento relativo dell’oscillare viaggiante (b), v(t) , e dell’accelerazione assoluta dell” oscillato- 
re viaggiante (c), ¥,(t)=g—@, v(1)- 26, @, W(t). Anche in questo caso i risultati forniti dal metodo 
MSE sono a svantaggio di sicurezza. 

La tecnica proposta si presta ad essere efficacemente applicata per valutare, attraverso simulazioni 
Monte Carlo, le statistiche della risposta di oscillatore e trave in presenza di irregolarità nel contatto. In 
una seconda fase, quindi, è stata introdotto nel sistema l’irregolarità aleatoria r(t) , descritta probabilisti- 
camente dallo spettro di equazione (19), caratterizzato da A. = 31.25 m° e y, =0.4 rad/m. In Figura 7a 
sono rappresentate le densità di probabilità dell’abbassamento della trave in corrispondenza 


Tabella 2. Frequenze circolari e rapporti di smorzamento modali 
modo proposto tecnica MSE 

1 0 130 rad/s Ó, 194 rad/s Č 0.0465 

2 i da 148" È, 205" È, 0.0423 

3 |} 050 23 " È 6 249 " | & 0.0326 

4 Oa 308" O, 340 " È, 0.0250 

5 i ws 455" ©; 477." È; 0.0218 


Figura 4. Forme modali della trave. 
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| — proposto 
4 --- MSE 


u(x,t) [m] 


x [m] x [m] x [m] 
Figura 5. Deformata della trave negli istanti t= 0.25 s (a), t= 0.50 s (b), t= 0.75 s (c). 


dell’oscillatore viaggiante, per x,(t)= 0.25 £, , 0.50L, e 0.75L,. Le densità di probabilità sono state 
costruite generando 1500 campioni del processo forzante N(t) attraverso la formula di equazione (17), 
ed applicando per ciascuno di essi lo schema numerico di equazione (26). La figura evidenzia che il coef- 
ficiente di variazione è elevato solo quando l’oscillatore viaggiante è in mezzeria. In Figura 7b sono con- 
frontate le densità di probabilità fornite dal metodo proposto e dal metodo MSE per l’abbassamento della 
trave quando x,(1)=0.50 L, : si osserva che il metodo MSE predice media e varianza inferiori a quelle 
del metodo proposto. In Figura 7c, infine, sono confrontate le densità di probabilità fornite dal metodo 
proposto e dal metodo MSE per il picco massimo dell’accelerazione assoluta dell’oscillatore viaggiante. 
In questo caso, per il sistema considerato nelle applicazioni numeriche, si rileva che i due metodi risultano 
in ottimo accordo, anche in corrispondenza delle code. 


7. CONCLUSIONI 
Nella memoria è stata presentata una tecnica che consente di eseguire l’analisi nel dominio del tempo del- 
le vibrazioni indotte da un oscillatore ad 1 Gdl che percorre una trave elastica, modellata al continuo, su 
un suolo viscoelastico. 


—— proposto 


-0.0001 95-77 


o 05 1 
t [s] 


0 0.5 1.5 


1 
t [s] 


1.5 


Figura 6. Storie temporali dell’abbassamento della trave in mezzeria (a), dell’abbassamento relativo 
dell’oscillatore viaggiante (b) e della sua accelerazione assoluta (c). 
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160000 50000 103 
| |12025L gy i MSE b) c) 
rg 120000- T y T b 
J I proposto 
= x,=0.75 L, = 30000] ʻi y 
= 800004 = fot 2, 014 
= | = 200004 14 ~ 3 
o -oso| £ i) dii 
E, 40000, UU La, & 10000 1 ) E 0.013 — proposto 
7 A ZA Sie MSE 
oth ot 0.001 1—— 
0 0.0002 0.0004 0.0006 0.0002 0.0004 0.0006 0.2 0.6 1 14 
u [m] u [m] y, [m/s] 


Figura 7. Densità di probabilità dell’abbassamento della trave (a e b) e del picco massimo dell’accelerazione 
dell’oscillatore viaggiante (c). 


Rispetto alle formulazioni disponili in letteratura, in cui si utilizzano valori equivalenti di rigidezza e 
dissipazione della fondazione, la tecnica proposta consente di studiare l’interazione oscillatore-trave adot- 
tando modelli reologici più accurati per il suolo viscoelastico. In particolare, nella memoria si è fatto rife- 
rimento al cosiddetto modello di Zener, che può essere efficacemente impiegato per il materiale elastome- 
rico da inserire tra rotaia e sottopiastra per ridurre le vibrazioni indotte dal transito dei veicoli ferroviari. 
Per questo modello sono state presentate le equazioni di stato, che sono state poi accoppiate con le equa- 
zioni del moto dell’oscillatore viaggiante e della trave elastica. 

Nelle applicazioni numeriche è stata evidenziata l’entità dell’errore che si commette utilizzando valo- 
ri equivalenti per la rigidezza e la dissipazione dell’elastomero (modello di Kelvin-Voigt equivalente). E° 
stata inoltre verificata la possibilità di impiegare la tecnica proposta per valutare le statistiche della rispo- 
sta in presenza di un’irregolarità aleatoria nel contatto oscillatore-trave. Per il sistema considerato nelle 
applicazioni, in particolare, le statistiche del picco massimo dell’accelerazione assoluta dell’oscillatore 
viaggiante possono essere accuratamente valutate anche con il modello di Kelvin-Voigt equivalente, men- 
tre per le statistiche della risposta della trave è necessario utilizzare modelli più accurati. 
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